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SUMMARY 


The process of performing an automated aeroelastic stability analysis for an 
elastic-bladed helicopter rotor is discussed. A symbolic manipulation program, 
written in FORTRAN, is used to aid in the derivation of the governing equations of 
motion for the elastic-bladed rotor. The blades undergo coupled bending and tor- 
sional deformations. Two-dimensional quasi-steady aerodynamics below stall are 
used. Although reversed flow effects are neglected, unsteady effects, modeled as 
dynamic inflow are included. Using a Lagrangian approach, the governing equations 
are derived in generalized coordinates using the symbolic program. The symbolic 
program generates the steady and perturbed equations and writes into subroutines to 
be called by numerical routines. The symbolic program can operate both on expres- 
sions and on matrices. For the case of hovering flight, the blade and dynamic 
inflow equations are converted to equations in a multiblade coordinate system by 
rearranging the coefficients of the equations. For the case of forward flight, the 
multiblade equations are obtained through the symbolic program. The final multi- 
blade equations are capable of accommodating any number of elastic blade modes. The 
computer implementation of this procedure consists of three stages: (1) the sym- 

bolic derivation of equations; (2) the coding of the equations into subroutines; and 
(3) the numerical study after identifying mass, damping, and stiffness coefficients 
for each equation. Damping results are presented in hover and in forward flight 
with, and without, dynamic inflow effects for various rotor blade models, including 
rigid blade lag-flap, elastic flap-lag, flap-lag-torsion, and quasi-static tor- 
sion. Results from dynamic inflow effects which are obtained from a lift deficiency 
function for a quasi-static inflow model in hover are also presented. 

The numerical results for hovering flight indicate that dynamic inflow 
increases the lead-lag regress ing-mode damping for torsionally rigid blades. For 
torsionally flexible blades, the dynamic inflow effect depends on the elastic cou- 
pling parameter. For zero elastic coupling, dynamic inflow increases the modal 
damping. For full elastic coupling, it decreases the damping. This implies that 
there exists an elastic coupling parameter value for which dynamic inflow effects 
are negligible. The study for forward flight indicates that for a large number of 
degrees of freedom and nonlinear models, the amount of input data to the symbolic 
program increases exponentially, making it inconvenient to explicitly consider the 
harmonic and multiblade equations. However, a combination of symbolic and numerical 
programs at the proper stage in the derivation process makes it effective and bene- 
ficial to obtain the stability results from this approach. The numerical study 
indicates that dynamic inflow does change the magnitude of the predicted damping, 
yet its influence on damping trends is generally small with varying advance ratio or 
elastic-coupling parameter for torsionally flexible blades. In this report, Part I 
describes the symbolic program concepts. Part II presents the numerical results 
obtained using this process. 
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INTRODUCTION 


It is a general experience that derivation of governing equations of motion for 
complex structures represents a task of significant magnitude and is subject to 
error when performed by hand. When consideration is given to helicopter aeroelastic 
problems in particular, even for a rigid flap-lag model, the derivation of the 
governing equations of motion and the multiblade coordinate transformation requires 
much time and determination of accuracy by independent means. Experience with 
elastic hingeless rotor blade analyses indicates that for a given ordering scheme, 
the final equations differ in small nonlinear terms in the process of derivation 
depending on the stage at which the ordering scheme is applied. When the ordering 
scheme is consistently applied at a later stage in the derivation process, in gen- 
eral, this process requires more time and more lengthy independent checking. 

In this situation, a symbolic manipulator allows the analyst to share the 
algebra with the computer. Use of symbolic programs has been reported in several 
branches of science and engineering for over 30 years. There are many general 
symbolic processors available, for example, MACSYMA (ref. 1), REDUCE, etc., written 
in LISP. MACSYMA has been used in many structural applications (ref. 2), and REDUCE 
has been used in helicopter applications (ref. 3). When applied to a particular 
problem, the numerous capabilities available from such general purpose programs tend 
to slow the program execution and puts restrictions on computer memory. 

To circumvent these shortcomings, programs have been developed which are spe- 
cially tailored to given tasks, as for example, in celestial mechanics 
(refs. 4-9). These programs take advantage of the special form of the expressions 
to be manipulated. They deal with the powers of the fixed number of variables 
forming the expressions and they manipulate the algebraic operations. These pro- 
grams were written partly in FORTRAN and partly in machine language. In refer- 
ence 10, concepts have been presented to manipulate matrices and series (expres- 
sions) of general form using FORTRAN. Using these basic concepts, a program called 
HESL has been written in FORTRAN (ref. 11) and was used in deriving rigid blade 
helicopter rotor equations for aeroelastic analysis. 

In the present study, HESL has been extended and used in the derivation of 
elastic blade equations in generalized coordinates using a Lagrangian formulation. 
The equations can be derived for any given ordering scheme. These equations are 
coded into FORTRAN subroutines. The statements such as the subroutine names COMMON, 
DIMENSION, etc., which are required in making these subroutines are read as data 
during the coding. The program writes two subroutines, one for nonlinear equations 
to calculate the rotor trim position and the other for use with linearized perturba- 
tion equations to analyze stability. These subroutines are subsequently called by 
numerical routines which identify the mass, damping, stiffness, and load terms for 
each equation to form the required matrices which govern the blade behavior. Har- 
monic balance equations and transformed multiblade coordinate equations are obtained 
from the symbolic program. A schematic of the entire operation is shown in 
figure 1(a). The figure shows three segments of the analysis process; segment 1 for 
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computer derivation of the equations, segment 2 for incorporating the matrix ele- 
ments into FORTRAN subroutines, and segment 3 for numerical study. 

In this report, Part I describes the symbolic program aspects and Part II 
presents numerical results on the aeroelastic stability of an elastic rotor with 
dynamic inflow both in hover and in forward flight. 


I: SYMBOLIC PROCESSOR PROGRAM PRINCIPLES 


In this section, the concepts used to develop a symbolic (algebraic) manipula- 
tion program using FORTRAN are described. The basic manipulation of multiplying two 
expressions symbolically is then explained. Using this manipulation, the algebraic 
operations required in the aeroelastic stability and in response analysis of heli- 
copter rotors are developed and presented. The matrix operations are performed by 
assuming that the elements of the matrix are expressions. The remaining sections 
explain the mode of input to the symbolic program. A sample output from the compu- 
ter program is given in appendix A. All the symbols used in this section are given 
in table A1 . 


Program HESL 

Reading and storing i nformation - The basic principle (refs. 10 and 11) is to 
associate numbers with the variables (blade deflections, rotations, and time deriva- 
tives) that are to be manipulated. This allows numbers instead of variables to be 
manipulated. The program automatically assigns a number to the variable as soon as 
it reads it for the first time. Initially the program assumes that an expression or 
a relation is composed of a number of individual terms. Each term within an expres- 
sion or a relation consists of (1) a single numerical coefficient, and (2) a pattern 
which consists of the numbers associated with each variable in the term. For exam- 
ple, tension T is given by 


T = EA u' + £ v' 2 + | w' 2 ( la) 

where E is Young's modulus, A is the area of cross section, u' is the derivative 
of the extensional displacement, and v' and w' are bending slopes in two planes. 
The equation for T could be defined as an expression given by 

T = E*A*( 1 . 0*US + 0.5*VS*VS + 0.5*WS*WS) (1b) 

(expression) = (term) + (term) + (term) 

which contains three terms. Each term has a numerical coefficient and is made up of 

variables US, VS, WS, E, and A. The data for each term are read in alphanumeric 

format (e.g., F10.0, A1, A4), and numbers are assigned by the program to variables 
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E, A, US, VS, and WS. Let 1, 2, 3, 4, and 5 respectively, be the numbers assigned 
to these variables. Then the first term (E A US) has a numerical coefficient of 1.0 
and a pattern of 1 2 3. The second term (0.5 E A VS VS) has a numerical coefficient 
of 0.5 and a pattern of 1 2 4 4; the third term has a numerical coefficient of 0.5 
and a pattern of 1 2 5 5. In the program, the expression T is (where brackets are 
used for convenience) 

T = 1 ,0( 1 )(2)(3) + 0.5(1)(2)(4)(4) + 0.5(1 )(2)(5) (5) (2) 

As another example, a relation given in a table which is used for subsequent substi- 
tution later in the derivation process might be 

9 r 9 + 9 cos i|> + 0 sin ip 

0 c s 

THTA = 1 . 0*THT0 + 1 .0*THTC*CSCY + 1 .0»THTS*SNCY (3) 

(relation) = (term) + (term) + (term) 

which shows the total pitch angle as given by collective and cyclic pitch angles 
and ip is nondimens ionalized time. By letting 6, 7, 8, 9, 10, 11 represent THETA, 
THT0, THIC, CSCY, THTS, and SNCY, equation (3) in the program is 

(6) = 1.0(7) + 1 .0(8) (9) + 1.0(10)01) (4) 

The numerical coefficient and the pattern of each term are stored in a coefficient- 
storage stack (or array), and in a pattern storage stack. These stacks are common 
to all expressions and to all relations given for integration, differentiation, 
perturbation, trigonometric relations for multiblade coordinate transformation, 
etc. The definition of expressions and relations and their storing is schematically 
shown in figure 1(b). A separate array (expression array) is provided to store the 
number of terms of the expression, such as T, its identification number, and the 
starting position in both the storage stacks. This information is required for 
subsequent expression retrieval. Similarly, the relations (right-hand side of THTA, 
in eq. 3) required for integration, differentiation, perturbation, etc., are read as 
tables of relations and are stored in stacks. The starting position and number of 
terms of the relations (THTA) for all the tables (distinguished by different table 
names), table identification number and left-hand side of the relation (THTA), are 
stored in a table array for retrieval. If the relations are not available in closed 
form, numerical schemes are used. The elements of a matrix are assumed to be formed 
by expressions. A matrix array stores the number of rows and columns of the 
matrices and the position of the corresponding numerical coefficient and pattern of 
the elements of the matrix in the storage stacks. 

Special symbols are provided to identify a variable, symbol ' ' (blank), an 
expression, symbol , a matrix, symbol a table containing known relations, 

symbol For example, expression T (eq. (1)) is identified as an expression in 

the program by placing the symbol % in front of the expression. Similarly, varia- 
bles WS, US, VS, THT0, THTC, THTS, CSCY, and SNCY are recognized as variables by 
placing a blank before them. Let the table containing relation (eq. (3)) be named 
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TRIG. It is identified in the program by placing '©' in front of the table name as 
©TRIG. 


Algebraic manipulation - The basic manipulation is the multiplication of two 
expressions, say X and Y, to get expression Z. The details of X and Y are 
obtained from the expression array (number of terms, starting location in common 
stacks), and from the common stacks (numerical coefficient and pattern of each 
term). The multiplication is then performed term by term. The product of two terms 
yields a new term. The numerical coefficient of the new term is the product of the 
numerical coefficients of the two terms. The pattern of the new term is the conca- 
tenation of the patterns of the two terms in the product. For example, if term 
E A US and term 0.5E A VS VS are multiplied, we get the new term 
0.5 E A US E A VS VS, which in symbolic notation becomes 

0.5(1)(2)(3)(1)(2)(4)(l|) 

After all the terms are multiplied, the resulting terms are stored in the 
common stacks. The location of these terms in the common stacks, and the number of 
terms of the expression Z, and identification number for Z are stored in the 
expression array for retrieval. 

The second manipulation is substituting relations, given in a table within the 
computer code, into an expression. This is conceived as substituting the pattern 
for the relation (such as the pattern for THTA) into an expression containing THTA 
and carrying out the required multiplication, as explained above. The position of 
THTA in the common stacks is obtained from the table array. It can be seen that 
integration, perturbation, multiblade coordinate transformation, etc., can be easily 
performed with the known relations in this way. For differentiation, the same 
substitution technique is used; however, care is taken to store the dependent and 
independent variables separately. It can be realized that substitution of known 
relations for integration and differentiation is similar to the table look up proce- 
dure used in many aerodynamic calculations. 

A compacting capability is provided for the addition of similar terms within an 
expression. This is done by comparing two terms within an expression for patterns, 
and, if they are identical, numerical coefficients are added and the number of terms 
in the expression is reduced by one. 

The ordering scheme is applied in the following manner. The order of each 
variable and the overall ordering scheme are read as user-supplied data. Once all 
the terms in each expression are determined, the order of each variable is added on 
a term-by-term basis. If the order sum for a particular term is greater than the 
specified order to be retained, the term is neglected. 

In order to collect all the coefficients that multiply a state variable in the 
final equations, the pattern of each term is checked to see whether the number 
assigned to that state variable is present. If it is present, the term is stored 
separately . 
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To perform a multiblade transformation, the relevant trigonometric relations 
are supplied as a table of relations. The multiblade functions 1, cos *|>, sin *|>, 
etc., are read as data depending on the number of blades. For example, these func- 
tions for a f ive-bladed rotor are 1, cos 4>, sin 41 , cos 2*1), and sin 2*1). The multi- 
blade expansion of each generalized coordinate is given as data. Then the required 
multiplications and substitutions are performed. The terms are then checked for 
cos N*|) and sin N*ii (N = number of blades). Terms which do not have integer multiples 
of number of blades are deleted. Harmonic balance equations are obtained by substi- 
tuting the trigonometric relations and by collecting all like harmonics. 

Using the above concepts, the program HESL was written in FORTRAN and used to 
generate the governing equations of motion of a coupled rigid fuselage and rigid 
rotor blade analysis (ref. 11). HESL can perform all the algebraic manipulations 
required in the equation derivation. It can perform these operations at an expres- 
sion level or matrix level. Matrix multiplication, transpose multiplication, addi- 
tion, subtraction, and multiplication of matrix elements with a constant or expres- 
sion can be done by using HESL. The matrix operations are done by assuming that the 
elements of the matrix are expressions. The matrix elements can be changed to 
expressions. As explained above, integration, differentiation, perturbation, and 
multiblade coordinate transformations are done by substituting known relations 
provided as data to the program. Numerical schemes are used if necessary. It 
should be noted that for the harmonic balance equations and multiblade coordinate 
transformation, trigonometric relations (products of sines and cosines expressed as 
sums of sines and cosines) are given as a table of relations and substitution is 
done. Equations can be derived for any ordering scheme. 

The program has 40 subroutines excluding numerical routines, and are called 
through a main program by reading the manipulations to be performed as commands. 

For example, the command FORM EXPRESSION is used to call the necessary subroutines 
to multiply two or more expressions. The command FORM MATRIX is used to call the 
subroutines to multiply two matrices. The command SUBSTITUTE TABLE INTO EXPRESSION 
will substitute known relations into an expression, etc. The HESL code has 
22 recognized commands. Table 1 lists each of the various commands and their 
functions. 

The basic input to the program are the position vector, transformation rela- 
tions, order of each variable, and the ordering scheme. If the relations are known, 
the relations for integration, differentiation, perturbation, etc., are given as 
tables of relations. The nonlinear equations are obtained by calculating strain 
energy, kinetic energy, and work done in generalized coordinates, then a Lagrangian 
formulation is used to obtain the governing equations of motion. The equations are 
linearized using the perturbation relations. When using trigonometric relations, 
harmonic balance equations and multiblade equations are obtained. The multiblade 
summation rules are embedded in the program. It should be noted that apart from the 
basic relations, the data to generate a required function has to be given as 
input. For example, one has to define what functions are to be multiplied to get 
strain energy. For a different blade model only the basic relations that define the 
change in the model need to be defined. 
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Data and program aspects - The main program initializes the program data, reads 
the commands, and calls the appropriate subroutines. When it encounters the com- 
mand, END OF DATA, execution is stopped. In order to keep the available core space 
to a maximum, all the expressions which are no longer needed are erased from the 
common stacks by a RESET COUNTER. 

The first input line in the program deck defines the identifiers to be used for 
a blank variable 'bbb', a variable 1 b 1 , an expression * % ' , a table '§', and a matrix 
in that order. 

This first input line is immediately followed by the input for the algebraic 
manipulations : 


12345678 column numbers (not an input) 
bbbbb£@# (b denotes a blank space) 

In general the input is the command name, followed by the name of the expression 
(table, matrix) and number of terms (relations, size) of the expression (table, 
matrix) followed by the term (relation, term) details. These expressions are in 
fixed format so as to provide consistency and to avoid confusion. 

The names of variables, expressions, matrices, and tables are restricted to 
four alphanumeric characters and are read in A4 format. The identifiers are read in 
A1 format. For example, variables WS and VS are recognized by reading them as bbbWS 
and bbbVS. Expressions T and VEL are read as ?bbbT and ?bVEL, matrix LAFP is read 
as #LAFP. Table TRIG which contains given relations is read as @TRIG. Other sym- 
bols are used to separate them from this predefined set. For example, symbol in 
*E2D1 is used for identifying an ordering scheme, and symbol "e' in ePECF is used 
for identifying a group having certain variables. 

The terms of an expression or relation are read as input with only one term for 
each input line. The first ten columns of the line are reserved for the numerical 
constant and the remaining columns are used for the variables (or expressions) 
forming the term. The format is 


F10.0,14(A1,A4) (5) 

For example, the term 0.5 WS WS VS VS is read as 

123456789012345678901234567890 column numbers (not an input) 

0 . 5bbbbbbbbbbWSbbbWSbbbVSbbbVS 

The input and output of the program can be best explained for an example by 
deriving equations for the flap motion of a rotor blade model. Appendix A gives the 
complete input and output for the problem chosen. The problem definition starts 
with a flap-lag transformation definition, but is subsequently reduced to a flap- 
only model. This simplified model is used to clarify the program aspects rather 
than clarifying the modeling aspects. It should be noted that other forms of 
achieving the required objective may exist in addition to those defined here which 
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use the basic algebraic manipulations. Table A1 shows all the symbols used in the 
program. 

The following explains how commands and data are read to the program. The data 
can be considered as consisting of two parts. One is directly read by the read 
commands (READ EXPRESSION, READ TABLE, etc.). For the second, data operations for 
the manipulations, i.e., FORM EXPRESSION, FORM MATRIX, SUBSTITUTE TABLE INTO EXPRES- 
SION, etc., are used. The commands can be divided into three parts: (1) Input 

Commands; (b) Algebraic Manipulation Commands; and (3) Application Commands (see 
table 1). 


General Rules 

1. Commands names are entered from the first column. Only the first eight 
characters of the command are important for the operation. 

2. The names of variables, expressions, matrices, and tables are always asso- 
ciated with their identifiers, ' ', and '§', respectively (read in the 

first input line). 

3. Other identifiers can be used to distinguish them from variables, expres- 
sions, matrices, and tables, for example, *,e,1, to be used in specifying ordering 
schemes, group names for collection of coefficients. 

Input commands - These commands are used to read the input data, such as details 
of the terms forming the expressions, matrices, and tables of relations. 

Command "READ EXPRESSION". This command is used to input an expression con- 
sisting of variables only. The input sequence is: 

1. Command name; 

2. Expression identifier (%) , expression name, and the number of terms in the 
expression in format ( A 1 , A4 , 12 ) ; 

3. For each term, input details of the term in format (F10.0, 1 4 ( A 1 , A4) ) . 
Example to input the strain energy of the blade given by 


where K 0 is the flapping stiffness and 8 is the flapping angle and to read 
expression T given by 


T = EA u' 
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Input 


Explanation 


123456789012345678901234567890 

READ EXPRESSION 

5SSEBL01 

0.5 KBT BT BT 


READ EXPRESSION 
% T03 

1 .0 E 

0.5 E 

0.5 E 


A US 
A WS WS 

A VS VS 


(column numbers, not an input) 
command name. 

identifier, expression name, number of terms, 
numerical coeff . , variables KBT,BT forming the 
term. 

command name. 

identifier, expression name, number of terms, 
numerical coeffi., variables E,A, US. 
numerical coeffi., variables E,A, WS,WS. 
numerical coeffi., variables E,A, VS, VS. 


Command "READ MATRIX". This command is used to input the elements of a 
matrix. The term details of the matrix A are read as ( A( I , J) , J= 1 ,N) , 1= 1 ,M) . The 
terms of the matrix element can be variables or already defined expressions or a 
combination of these expressions. The input sequence is: 


1 . Command name ; 

2. Matrix identifier (#), matrix name, size of the matrix (rows by columns), 
and ordering scheme if required, in the format ( A 1 ,A4,2I2, A1 , A4) ; 

3. Then for each element of the matrix, read the number of terms and for each 
term read the term details in format (F10.0, 1 4 ( A 1 , A4) . Example: To input the flap 
transformation matrix, given as 

sin 8 
0 

cos B_ 

where 8 is the flapping angle. 

Input Explanation 


flp 


’cos 8 0 

0 1 
sin 8 0 


12345678901234567890 
1234567890 
READ MATRIX 
#TFLP0303 

01 

1.0 C0SB 
01 


(column numbers, not an input) 
command name. 

matrix identifier, matrix name, number of rows, 
number of columns, 
number of terms in element T( 1,1). 
numerical coefficient, variable C0SB 
(C0SB is cos 8). 

number of terms in element T(1,2). 
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Input (cont'd) 


Explanation 


0.0 


01 


-1.0 

SINB 

01 


0.0 


01 


1.0 


01 


0.0 


01 


1.0 

SINB 

01 


0.0 


01 


1 .0 

COSB 


Command "1 


numerical coefficient, 
number of terms in element T(1,3). 
numerical coefficient, variable SINB 
(SINB is sin s). 

number of terms in element T(2,1). 


D TABLE FOR SUBSTITUTION". This command is used for reading a 
table of relations which is used in subsequent substitutions. Each table of rela- 
tions is identified by a table name. HESL identifies two types of tables of rela- 
tions. Type 1 is a substitution of relations of the type v = v + 6v, and type 2 is 
a relations table which has powers (>1) and a product of the variables such as 


sin* 1 ip = 0.5-O.5*cos 2ip. 
form of 
their powers. 


A S, B m C n = 


In type 2, the relations can be handled in the general 


terms, where A, B, , 
For example, sin" 5 ip cos £ 


nd C are variables, and l, m, and n are 
e sin^ $ = terms. The input sequence is: 


1 . Command name ; 

2. Table identifier (#), table name, number of relations in the table, and 
type of the table used in format (A1,A4,I2). A zero value indicates type 1 rela-- 
tions, a nonzero value indicates type 2 relations; 

3. Then, for each relation, input: 

(a) the number of terms on the right-hand side of the relation, the left- 
hand side variables and their power in format ( 12 ,4( A4, 12) ) ; 

(b) the details of the terms (of the right-hand side of the relation) in 
format (F10.0, 14( A1 , A4) ) . 

Example: To input two tables of relations. Table 1, named SUPR contains 

powers of unity of the variables in the relation. Table 2 contains powers greater 
than the unity of the variables and the product of variables. 


Table 1 (1) sin c = 0.0; (2) cos z, - 1.0; (3) 9 = 9 + 0 cos ip + 0 sin ip 

O C o 
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Table 2 

2 2 
(1) sin B = 1.0 - 1.0 cos 8 cos 8; (2) sin ip cos ip = 0.25 sin ip + 0.25 sin 3>P 

where r, is the lagging angle and 6 is the flapping angle, and ip is nondimen- 

sional time. 


Input 


Explanation 


123456789012345678901234567890 
READ TABLE FOR SUBSTITUTION 
@SUPR0300 

01 SINZ01 


0.0 

01 C0SZ01 
1.0 

03 THTA01 

1.0 THTO 

1.0 THTC CSCY 

1.0 THTS SNCY 

READ TABLE FOR SUBSTITUTION 

@DUMY020 1 

02 SINB02 
1 .0 

1.0 COSB COSB 
02 SNCY01 CSCY02 

0.25 SNCY 
0.25 S3CY 


(column numbers, not an input) 
command name. 

table identifier, table name, number of 
relations, table type. 

number of terms on the right hand side, relation 
name, its power. 

numerical coeff. (term details), 
same for second relation. 


command name. 

identifier, table name, number of terms, table 
type. 

number of terms on the right hand side, relation 
name, its power, 
numerical coefficient, 
numerical coefficient, variable COSB. 
number of terms, left hand side variables with 
their powers. 

numerical coefficient, variable SNCY. 
numerical coefficient, variable S3CY. 


Command "READ DIFFERENTIATION TABLE". This command is used for reading the 
rules of differentiation of variables. The input sequence is: 


1 . Command name ; 

2. Table identifier, differentiation relations table name and the number of 
relations in the table in format (A1,A4,I2). 


Then for each relation provide: 

(a) The independent variable for differentiation, the dependent variable to 
be differentiated, and the number of terms in the differentiation relation in format 
(2{A1,A4)I2); 

(b) Each term detail in format (F10.0, 1 4 ( A 1 , A4) ) ; 


1 1 



Example: To read a table of relations containing the differentiation rules, 

given by 


M _ i n - M _ i o- 3 cos g 
38 ' ' ’ gg ' ’ 38 


_ 3 sin 6 „ 3 8 

■ Sln 8 ! “ 5 s — = cos 6 ; 57 = a 


|S = S_S 22 J = . sln s .j; = cos B- 8 ; = -sin *; §7 sin * = cos* 


3x 


3t 


where 6 is the flapping angle, and denote time. 


Input 


Explanation 


1234567890123M567890 1234567890 
READ DIFFERENTIATION TABLE 
gDERVIO 

BT BT01 

1.0 

BTD BTD01 

1.0 

BT C0SB01 
-1.0 SINB 

BT SINB01 

1.0 COSB 
TAU BT01 

1.0 BTD 
TAU BTD01 

1.0 BTDD 
TAU C0SB01 
-1.0 SINB BTD 
TAU SINB01 

1.0 COSB BTD 
TAU CSCY01 
-1.0 SNCY 
TAU SNCY01 

1.0 CSCY 


(column numbers, not an input) 
command name. 

identifier, name of the table, number of 
relations. 

independent variable, dependent variable, its 
power. 

numerical coefficient. 


Command "READ GROUP AND ORDER OF THE VARIABLES". This command is used to input 
the variables group to which they belong and their order of magnitude. Assigning 
group numbers and order of magnitude to the variables shows their importance in the 
analysis. For this example, original variables and steady quantities, 

8 , 8 , 8 , C^ /a are considered to belong to group 1 and perturbed quantities 

68 , 58, 68 °belong to group 2. This allows retention of only linear terms in per- 
turbation quantities. The input sequence is: 
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1 . Command name ; 


2. The total number of variables to which a group and order of magnitude is 
assigned (highest order implies lowest importance) in format 13; 

3. For each variable, provide the identifier for variable, variable name, a 
blank and its group number and order of magnitude in format ( A 1 , A4, IX, 212) . Eight 
variables can be typed per line. 

Example: To input the following variables showing their group and order of 

magnitude (gr = group, 0 = order). 

6 gr ( 1 ) 0(1); 8 gr( 1 ) 0(1); 8 gr(1) 0(1); 

68 gr(2) 0(1); 68 gr(1) 0(1); 68 gr(2) 0(1); 

C d /a gr( 1 ) 0(2); 9 gr(1) 0(1) 


Input Explanation 

READ GROUP AND ORDER OF command name. 

THE VARIABLES 

b08 total number of variables assigned group and 

order . 

BTDD 0101 BTD 0101 BT 0101 DBDD 0201 DBD 0201 DB 0201 CD0A 0102 THTA 010 

Variables BTDD, BTD, BT are assigned group 1(01) and order 1(01), variable CD0A is 
assigned group 1(01) and order 2(02), variables DBDD, DBD, DB are assigned 
group 2(02) and order 1(01), since they are considered as perturbed quantities. 

Note that all the variables have their identifier (a blank) in front of them. 

Note: If the order and group of any variable is not defined, the variable is 

assigned to have an order of zero (highest importance) and belong to group 1. 

Command "READ ORDERING SCHEME". This command is used for specifying the order- 
ing scheme. It specifies up to what order group (1) variables and group (2) varia- 
bles should be retained. The input is: 

1 . Command name ; 

2. Identifier for ordering scheme (*), the ordering scheme name and number of 
groups of variables considered in format ( A 1 , A4 , 12 ) ; 

3. For each group, provide in sequence the group number and the highest order 
of magnitude to be retained, in format (412). 

Example: To specify two groups of variables and retain group 1 variables up to 
total order of magnitude 2 and group 2 variables up to total order of magnitude 1 


13 


Retain 


Gr(1) ♦ 0(2) and Gr(2) ♦ 0(1) in a term. 


Input Explanation 

READ ORDERING SCHEME command name. 

*E2D102 identifier for ordering scheme, name of the 

ordering scheme, number of groups. 

01020201 group 1(01), highest order to be retained (02), 

group 2(02), highest order to be retained (01). 

Command "READ VARIABLES FOR COLLECTION OF COEFFICIENTS". This command identi- 
fies the name of a group of variables for which the terms containing these variables 
have to be collected. The input is: 

1. Command name; 

2. Group identifier (e), the name of the group and number of variables in 
format (A1,A4,I2); 

3. The variables in format (16(A1,A4)). 

Example: (1) To input a group name containing variables 8, 8, 8, 9 for 

collecting terms containing these variables. subsequently. (2) To input a group 
name containing variables 9 q , 0^, 8 Q 
lecting terms containing these variables 

Input Comment 

READ VARIABLES FOR COLLECTION command name. 

OF COEFFICIENTS 

ePECF04 group identifier, its name, number of variables in 

the in the group. 

BTDD BTD BT THTA variables. 

READ VARIABLES FOR COLLECTION 
OF COEFFICIENTS 
eMUCF12 

TH0 TH1 TH2 BDD0 BDD1 BDD2 BD0 BD1 BD2 B0 B1 B2 

Note that the terms containing these variables will be subsequently collected. Also 
note that variables are typed along with their identifier (a blank). 

Command "READ TAPE". This command is used to read and separate steady and 
perturbed quantities and to write on the disk to be used by the symbolic program. 

The input is: 


8 . 8 


8 _ , 8 . 


’ w 1. u 2’ o’ 
subsequently. 


V B 1 


i or col- 
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1 . Command name ; 


2. Expression name to be read from tape with its identifier, the number of the 
tape (disk) from which data has to be read; a nonzero value if perturbed and steady 
terms have to be separated, number of perturbation terms, number of the disk on 
which steady and perturbed terms have to be written, number of equations to be read 
in format ( A 1 , A4, 1012) . 

Example: To write the details of expression LAEQ on tape number 52. 

Input Explanation 

READ TAPE command name. 

56LAEQ5200 1 000000 1 expression identifier, expression name, number of 

the tape and other options. 

Algebraic manipulation commands - These commands are used to perform a single 
algebraic manipulation such as multiplication or substitution of relations for 
expressions and matrices. These commands mostly operate on the data read by input 
commands . 

Command "FORM EXPRESSION". This command is used to form another expression by 
addition, subtraction, and multiplication of variables and previously read expres- 
sions (by READ EXPRESSION command). The input is similar to that for READ EXPRES- 
SION command. The input sequence is: 

1 . Command name ; 

2. Expression identifier {%) , the name of the new expression to be formed, the 
number of terms forming the new expression, and if required, the name of the order- 
ing scheme, the name of the group containing the variables for which the coeffi- 
cients have to be collected, and a nonzero value to write on the disk (disk number), 
in format (A1,A4,I2,2(A1,A4),2I2); 

3. The details of each term in format (F10.0, 14( A1 ,A4) ) . 

Example: (1) To form a new expression c, from three already defined or read 
expressions, 5 ,H , with ordering scheme E2D1, and collect the terms containing 

variables in group^PE^F, to write on tape number 8 . (2) To form a new expression 
M g , with variables RAC2, e, r and already defined expressions U T and U p . 

Input 

FORM EXPRESSION command name. 

% RXD02*E2D1ePECFQ8 expression identifier, new expression name, number 

of terms, name of the ordering scheme, name of 
the group containing variables for collection, 
number to write on tape. 
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Input (cont'd) 

numerical coefficient, expression name, 
numerical coefficient, expressions RY and OMEZ. 


1.0 % RXT 

1.0 % RYjJOMEZ 


FORM EXPRESSION command name 

% MB02 expression MB with its identifier, number of terms 

-0.5 % UTJt UT THTA RAC2 RB numerical coefficient, expression UT, variables 

THTA, RAC2, RB. 

0.5 % UT% UP RAC2 RB similarly second term details. 

Command "FORM MATRIX". This command is used to multiply two matrices 
(C = A*B), including transpose multiplication (C = A transpose * B). The input 
sequence is: 

1 . Command name ; 

2. Names of the two matrices to be multiplied, (A,B), the name of the result- 
ing matrix (C) with matrix identifiers, name of the ordering scheme, and a nonzero 
value for transpose multiplication in format ( 4 ( A 1 , All ) 1 2 ) . 


Example: To multiply two matrices Tj a g and T^p to give T^g w: 

ordering scheme specified by E2D1 (if required transpose multiplication) 

[T flg' * tW IT flpl or ‘’’fig 1 = [T lag lT[T flp> 


Explanation 


FORM MATRIX command name. 

#TLAG#TFLP#TFLG*E2D100 matrix identifier, matrix name, matrix identifier, 

the resulting matrix name, matrix identifier, the 
resulting matrix name, ordering scheme, nonzero 
for transpose multiplication. 


Command "ADD MATRICES". This command is used to add or subtract two 
matrices. The input sequence is: 

1 . Command name ; 

2. Names of two matrices to be added or subtracted and the name of the 
resulting matrix with their identifiers and a zero for addition and nonzero for 
subtraction, in format ( 3 ( A 1 ,A4) ,12) . 

Example: To add matrices A and B to give matrix C 


[C] = [ A]±[B] 
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Input 


Explanation 


ADD MATRICES command name. 

# A# B# COO matrix identifier, matrix name, matrix identifier, 

matrix name, matrix identifier, matrix name, zero 
for addition (nonzero for subtraction). 

Command "CONSTANT MATRIX". This command is used for multiplying a matrix with 
a constant or an expression having only one term (both of which are read earlier in 
the program) . 

After multiplication the resultant matrix is written in the same location. For 
multiplying with expressions having more terms, define a diagonal matrix with this 
expression and perform matrix multiplication. The constant has to be specified as 
an expression before multiplication. The input sequence is: 


1 . Command name ; 

2. The matrix name, and the expression name with their corresponding identi- 
fiers in format 2(A1,A4). 

Example: To perform [A] = const * [A]. 

Input Explanation 


CONSTANT MATRIX 
# ANCONS 


command name. 

matrix identifier, matrix name, expression 
identifier, expression name. 


Command "MATRIX EXPRESSION". This command is used to convert matrix elements 
(usually a column vector) to expressions. The input sequence is: 

1 . Command name ; 

2. The number of expressions to be formed, the name of the matrix with its 
identifier ; 

3. For each expression, name of the expression with its identifier and its 
location in the matrix. 

Example: To redefine U T as given by element (2,1) of matrix VX2Y , and U p 

as element (3,1) of matrix VX2Y. 


Input Explanation 

MATRIX TO EXPRESSION command name. 

02#VX2Y number of expressions, matrix identifier and matrix 

name . 
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Input (cont'd) 


Explanation 


% UT02015S UP0301 expression identifier, name, element (2,1) and 

expression identifier, name, element (3,1). 

Command "SUBSTITUTE TABLE INTO EXPRESSION". This command is used to substitute 
relations given in a table (read earlier by READ TABLE FOR SUBSTITUTION command) 
into an expression with, or without, an ordering scheme. After substitution, coef- 
ficients for certain variables can be collected, and if required, can be written on 
a disk. The input sequence is: 

1 . Command name ; 

2. Name of the table containing the relations, the name of the expression, 
name of the resulting expression, ordering scheme, group for collecting coefficients 
with their respective identifiers, a nonzero value to write on the disk, in format 
(5(A1,A4),I2). 

Example: To substitute relations given in table SUPR in expression R1, result- 
ing in expression RX with an ordering scheme E2D1, group name for collecting PECF 
and no writing on the tape. 

Input Explanation 

SUBSTITUTE TABLE INTO command name. 

EXPRESSION 

§SUPR/6 RlyJ RX*E2D1 ePECFOO table name, old expression R1, resulting expression 

RX, ordering scheme E2D1, group for collection 
PECF with their respective identifiers, no 
writing on disk. 

Note: If E2D1 and PECF are not used, then leave blank spaces. 

Command "DIFFERENTIATE EXPRESSION". This command is used to differentiate an 
expression with respect to a specified variable using the differentiation rules 
previously read in the READ DIFFERENTIATION TABLE command (e.g., table DERV). The 
input sequence is: 

1 . Command name ; 

2. The variable name with which the expression has to be differentiated, the 
name of the expression to be differentiated, and the name of the resulting expres- 
sion in format ( 3 ( A 1 , A4 ) ) . 

Example: To differentiate expression r x (RX) with respect to item (TAU) result- 

ing in expression r xt (RXT). 
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Input 


Explanation 


DIFFERENTIATE EXPRESSION command name. 

TAUjt RX? RXT variable identifier, dependent variable, original 

expression and resulting expression with their 
identifiers. 

Application commands - Application commands are those which combine all, or some 
of the functions of the manipulation commands defined above to meet specific deriva- 
tion requirements. Application of Lagrangian, multiblade coordinate transformation, 
etc., require a number of operations to be done simultaneously. These commands are 
written to satisfy these requirements. 

Command "FORM LAGRANGIAN". By a careful consideration of the governing equa- 
tion deviation process, the required algebraic process can be suitably divided into 
a number of computational steps to be calculated without intermediate expression 
swell. For example if I is given by 

1 = J It [A * B1 * C dr (6a) 


It can be divided into 


I 


/ It * c , 


dr + 


/ * c 2 dr 


(6b) 


where A, B, and C, etc., are expressions and t is the independent variable. Here 
the evaluation of I is performed in two computational steps. The contribution 
from each computational step to the final expression is added to form the final 
expression. Each computational step requires a definition of expressions A^, B^, 
and C^. It also requires that there be relevant integral tables for integration, a 
differentiation table for differentiation, an independent variable, and perturbation 
relations for linearization. The evaluation of the total expression may also 
require an ordering scheme, or a group name for collecting terms. Accordingly, the 
input consists of: (1) the number of computational steps into which the problem is 

divided, the maximum number of tables of relations used in the evaluation, identi- 
fiers (1 or 0) for suppressing terms according to 1 + = 1 are required, print- 

ing before collection of terms, and for writing FORTRAN subroutines on disk; (2) the 
name of the resulting expression from this command, the name of the ordering scheme, 
and the name of the group having variables for which the coefficients have to be 
collected; (3) then for each computational step provide the names of the three 
expressions involved in the computational step, the name of the variable with which 
the product of the first and second expression is to be differentiated (leave blank 
if differentiation is not needed), and the names of integration and perturbation 
tables. 
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Example: To evaluate 


beq 


[ 3 - 

J 3' 


[r r ] 
L x xq J 


c dx + 


/ 


cM c dx 

P 


This is equal to two computational steps. The tables of required relations are INRL 
for integration, PERT for perturbation, DYNM for integration in aerodynamic terms, 
the independent variable x (TAU), the final equation name, S b (SBEQ), the ordering 
scheme, E2D1 , and the group name for the collection of terms, PECF. 


Input 


Explanation 


FORM LAGRANGIAN 
0202000100 

ItSBEQe PECF*E2D 1 

% RXD/tRXDQ^CONS TAU§INRL@PERT 

$cons* mb /Icons @dynm@pert 


command name. 

two computational steps, using two tables of 
relations and other options, 
final expression name, group for collection, 
ordering scheme, 
expressions, variable, tables 
expressions, variables, tables 


Note: When comparing the above inputs to equation (6b), A is RXD, B is RXDQ, C 
is CONS, x is TAU, with integration, perturbation relation in INRL and PERT. A 
similar process is repeated for the second term. 

Command "INITIALIZE MULTIBLADE". In performing the multiblade transformation, 
the single-blade equation has to be converted from a rotating frame of reference to 
a nonrotating frame of reference. This requires that the single-blade equations be 
multiplied and multiblade functions 1, cos 4>> sin 4>, (-1), etc., depending on the 
number of blades. This command is used to read a number of blades and to generate 
these expressions. This has to be called before performing a multiblade coordinate 
transformation. The input is: 


1 . Command name ; 


2. The number of blades, and (-1), cos 4>, sin 4>> cos^ 4), 
sin^4», cos 1 * 4>, sin 1 * 4), cos^ 4>, sin^ ip, cos° 4>, sin° 4), cos' 4>, 
(12, 3X, 1 5 ( A 1 ,A4)) . 


. 2 3 

sin 4>, cos 4>> 

sin^ 4>, in format 


Example: To read multiblade functions for a three-bladed rotor also cos(34<) 

and sin^ 4». 


Input 


Comments 


INITIALIZE MULTIBLADE 
03 M0PK CSCY SNCY C2CY 
S2SY C3CY S3CY 


command name. 

number of blades, multiblade functions. 

cos(34>), sin(34<) • 
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Command "PERFORM MULTIBLADE TRANSFORMATION". This command is used to obtain 
equations in multiblade coordinates from the single blade equations. The input 
sequence is: 

1 . Command name ; 

2. Number of resulting multiblade equations (equal to number of blades). 

Note: For dynamic inflow equations, the number of equations is one; 

3. The names of the tables containing multiblade relations, trigonometric 
relations, the name of the group containing multiblade variables for collection of 
coefficients; 

4. The name of the expression to be transformed, the names of the resulting 
multiblade expressions. 

Example: To convert equation S^gq into multiblade equations for a three- 
bladed rotor. Multiblade expansion of each degree of freedom is given in table 
MULB, trigonometric relations are given in table TRIG, the variables for collection 
of coefficients is given in MUCF. The final multiblade equations are named AEM1, 
AEM2, and AEM3. 

Input 

PERFORM MULT I BLADE 
TRANSFORMATION 
03 

@MULB@TRIG&MUCF 
5JSBEQ* AEM1 *AEM2*AEM3 


Explanation 
command name. 

number of resulting multiblade equations, 
relation tables (MULB, TRIG), group name for 
collection (MUCF). 

single blade expression name and resulting 
multiblade expressions names with their 
identifiers. 


Command "WRITE TAPE". This command is used to write on disk or tape, the 
steady and perturbed equations. The input is: 

1 . Command name ; 

2. Name of the expression to be written on disk, and the number of the disk in 
format (A1,A4,I2). 

Example: To write expression ABCD on disk number 3. 

Input Explanation 

WRITE TAPE command name. 

%ABCD03 expression name, disk number. 
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Command "NEGATE EXPRESSION". This command is used to get the negative of an 
expression when required. The input is command name and expression name. 

Example: To rewrite (UT) = -lip 

Input Explanation 


NEGATE EXPRESSION command name. 

% UT expression name with identifier. 

Command "RESET COUNTER". This command is used to reset the counter (starting 
number of expressions, see reading and storing section), and to save specified 
expressions that are required for further manipulation and calculations, so that 
further computations can over-write the unneeded expressions and thus save memory. 
The input is: 

1 . Command name ; 

2. The name of the expression whose starting number will become the new start- 
ing number (see reading and storing), and number of expressions to be saved. 

3. The names of expressions to be saved. 

Example: To reset the counter from expression EXP5 and to save expressions 

EX25 and EX26, which are generated after EXP5 is generated. 


Input 


Explanation 


RESET COUNTER command name. 

?EXP502 identifier, name and number of expressions to be 

saved. 

£EX25ytEX26 identifier, expression name(s) to be saved. 

Command "END OF DATA". When this command is read, execution is stopped. This 
command can be placed at any location in the derivation process to check the program 
output. The input is: 

1 . Command name . 

Input Explanation 


END OF DATA 


command name. 


Appendix A shows the output for the equation derivation for a rigid rotor blade 
model. The problem initially formulated with flap-lag degrees of freedom is reduced 
to flap degree of freedom (by SUBSTITUTE TABLE INTO EXPRESSION command). The output 
includes both the single-blade equation SBEQ, and the complete multiblade equations, 
AEM1, AEM2, and AEM3 for a three-bladed rotor. During the equation generation, the 
required quantities are written on a separate disk (tape) area for subsequent 
FORTRAN translation. The program TRANS (see table 1) is executed to read this data 
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and code it into subroutines. The elements of each equation, and each coefficient 
of the state variables, defined in PECF or MUCF (for multiblade equations) are 
written as elements of a matrix. The coefficients are separated into mass, damping, 
and stiffness matrices by calling this subroutine through a numerical program. 

Comments 

An algebraic manipulation program using FORTRAN has been described. The pro- 
gram can operate on both expressions and on matrices. Because of the modularity of 
the program, other operations can be easily added. The following limitations are 
observed with regard to data input and the size of the individual terms. The pro- 
gram requires data as input for differentiation, integration (because of the table 
look-up procedure), and for other operations required in helicopter rotor problems 
such as perturbation and multiblade coordinate transformation. The amount of these 
data increase with the complexity of the problem. For example, the input data 
increase with the specified number of harmonics for a trim calculation, the number 
of blades, and the number of degrees of freedom describing the model. The program 
at present assumes that each term of an expression is made up of no more than 
20 variables. This assumption simplified the program development. If more varia- 
bles per term are expected, the dimensions need to be increased. 


II: NUMERICAL STUDY 


Introduction 

The hingeless rotor blade configuration reduces mechanical complexity and 
increases rotor control power and damping relative to rotors with articulated 
blades. However, it also introduces complex aeroelastic behavior that requires a 
rigorous analysis for effective design. The derivation of the equations governing 
this behavior is tedious and error-prone because of the nonlinear nature of the 
problem. This has led to the use of symbolic derivation and to automation of the 
entire process, from derivation to numerical calculation with only limited user 
interface required. In the present report, this automation of the aeroelastic 
analysis is applied to the flap, lag, and torsion dynamics of an elastic rotor blade 
in hover and in forward flight with dynamic inflow. 

A complete dynamic analysis of a hingeless rotor blade consists of all three 
elastic degrees of freedom--flap, lag, and torsion. Initial analyses focused on the 
investigation of flap-lag stability of torsionally rigid blades. The models con- 
sisted of a rigid blade with spring-restrained hinges at the hub to simulate bending 
flexibility. The stability of this type of model was analyzed for both hover 
(ref. 12) and forward flight (ref. 13). Flap-lag stability of elastic blades with 
uniform properties was studied by Ormiston and Hodges (ref. 12), based on a deriva- 
tion of nonlinear partial differential equations suitable for elastic hingeless 
blades. Similar equations were studied by Friedmann and Tong (ref. 14). 
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Concurrently with the flap-lag stability analyses, efforts were also made to 
investigate the complete blade problem by including blade torsional deflections. 
Friedmann and Tong (ref. 14) approximated the torsional deflection by rigid body 
pitching motion (root torsion); they found that torsion motion was important and 
that the stability characteristics were sensitive to the number and type of assumed 
bending mode shapes used. Flap-lag structural coupling was not included. Hodges 
and Ormiston (ref. 15) presented extensive numerical results for the stability 
characteristics of elastic hingeless blades with flap-lag-torsion motion in hover . 
They found that torsional deflections of hingeless rotor blades are strongly 
influenced by the nonlinear structural moments caused by flap and lead-lag 
bending. This bending-torsion structural coupling is proportional to the product of 
the flap and lead-lag bending curvatures and to the difference between the two 
bending flexibilities. This study also showed the effect of precone, structural 
coupling, and torsional rigidity on the isolated blade stability boundaries. 
Friedmann and Kottapalli (ref. 16) analyzed the coupled flap-lag-torsional dynamics 
of hingeless rotor blades in forward flight. They noted that nonlinearities are 
important in an aeroelastic stability analysis, and that forward flight is strongly 
coupled with the trim state. Only flapping motion was used in calculating the rotor 
trim condition. It was observed that forward flight is stabilizing for soft 
in-plane rotors and destabilizing for stiff in-plane rotors. In all these studies, 
the aerodynamic forces were obtained from strip theory based on a quasi-steady 
approximation of two-dimensional unsteady airfoil theory. 

Simultaneously, efforts have been made to improve the aerodynamic model used in 
these analyses by including unsteady airflow effects. One approach to include these 
effects is to model the induced velocity as a time-dependent, three degree-of- 
freedom system. This dynamic inflow theory has been applied to rigid blade flap-lag 
analyses both in hover and in forward flight (refs. 17-19), and to a coupled rotor- 
fuselage problem in hover (refs. 20 and 21). It was observed that the dynamic 
inflow increased the lag-regressing mode damping and reduced the body pitch and roll 
damping for the parameters considered. These analytical results correlated well 
with experimental results (ref. 20). However, the conclusions presented in refer- 
ences 17 to 21 were based on several restrictive assumptions: zero elastic cou- 

pling, fixed-solidity ratio, and rigid flap-lag blade models with no torsional 
flexibility. 

It was observed in reference 15 for hover and in reference 16 for forward 
flight that when compared with a coupled flap-lag-torsion analysis, flap-lag analy- 
ses underpredict the in-plane (lead-lag damping) damping. As a result, it was 
pointed out in reference 22 that for torsionally soft blades to correctly assess the 
influence of dynamic inflow, it is necessary to formulate a model with both elastic 
torsion dynamics and dynamic inflow effects. 

General nonlinear differential equations for an elastic rotor blade used in the 
above analyses have been developed by Hodges and Dowell (ref. 23), Kaza and 
Kvaternik (ref. 24), and Rosen and Friedmann (ref. 25). The models have elastic 
flap, lag, and torsion degrees of freedom. Nonlinearities owing to moderate deflec- 
tions are also included. In these studies it was observed that for a given ordering 


24 



scheme, the final equations differ in small nonlinear terms, depending on the stage 
in the derivation process at which the ordering scheme is applied. The application 
of the ordering scheme at a later stage in the derivation process requires greater 
time in deriving and in independently checking the final equations. This has led to 
attempts to share the algebra with computers through symbolic processors. With such 
a processor to derive the governing equations, it will also be possible to study the 
effect of different ordering schemes on the calculated results in a minimum amount 
of time. Both general purpose and special purpose programs have been developed and 
are available as mentioned in Part I of this report. The HESL program, which is 
appropriate to rotory-wing aeroelastic investigation, has been developed in 
reference 11. The principles and method of organization of the program are 
explained in Part I of this report. 

In the present report, the governing equations of motion are derived using 
HESL. The program generates the steady state and linearized perturbation equations 
in symbolic form and then codes them into FORTRAN subroutines. The coefficients for 
each equation and for each mode are identified through a numerical program. A 
Lagrangian formulation is used to obtain equations in generalized coordinates. The 
coupled flap-lag-torsion-dynamic inflow equations are converted to equations in a 
multiblade coordinate system by rearranging the coefficients in the case of hovering 
flight. The explicit multiblade equations are symbolically derived using HESL in 
the case of forward flight. The multiblade equations are capable of accommodating 
any number of elastic blade modes. The whole process, from derivation to numerical 
calculation, is automated with minimum user interface. 

The hingeless rotor aeroelastic stability results presented in this report 
using the symbolic program reflect the combined effect of an improved structural 
model (by including torsion) and an improved aerodynamic model (by including dynamic 
inflow). Data are presented for several analytical models, a rigid blade lag-flap 
analysis, an elastic blade flap-lag analysis, an elastic blade flap-lag-torsion 
analysis, and a modified flap-lag analysis in which torsional inertia and damping 
are neglected (quasi-torsion). Damping data with and without dynamic inflow and 
from a quasi-static inflow model, in which dynamic inflow effects can be included as 
a lift deficiency function without increasing the system dimension, are included. 

The following sections present the numerical study performed using both the 
symbolic program described in Part I and the numerical programs developed to analyze 
the resulting equations. The first part presents the procedure, solution methods, 
and numerical results obtained for hovering flight with and without dynamic 
inflow. The second part presents the same for forward flight. 

Hovering Flight 

Formulation - Figure 2 shows an elastic blade with the coordinate system used in 
this study. The blade has uniform mass and stiffness, no twist, and no chordwise 
offsets of the elastic axis, tension axis, or center of mass. The elastic axis is 
coincident with the x-axis of the x,y,z coordinate system which rotates with a 
constant angular velocity (P) about a fixed point at the origin. The y-axis lies 
in the plane of rotation, and the x-axis is rotated through a small angle (B nn ) 

pL. 
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from the plane of rotation. The deflections of the beam are u (axial), v (lagwise 
bending), and w (flapwise bending) of the elastic axis parallel to the x,y,z 
coordinates, respectively. A second coordinate system, x', y', and z' fixed to the 
blade, with y' and z' axes parallel to the beam cross section principal axes, 
moves with the blade as it undergoes bending displacements, torsional displacements, 
and pitch angle (0) rotation. Before deformation, the principal axes of the blade 
are rotated with respect to the undeformed coordinates by the pitch angle. After 
deformation, the elastic axis is displaced by u,v,w, and the blade is twisted 
through the angle <t>. The inflow dynamics couple with the blade dynamics as a 
feedback loop, as shown in figure 3. The total inflow (v^ is assumed to consist of 
a steady value (X), and dynamic inflow components (v n , v , and v ) which vary with 
time. 


In this study the entire problem formulation is performed by the computer with 
minimum user interface other than a specification of blade geometry and desired 
blade model representation. In general, the formulation of the rotary-wing aero- 
elastic problem consists of the following steps (see refs. 22-25 for greater 
detail): writing the transformation matrices between the coordinate systems before 

and after deformation; calculating the position vector of a mass point of the 
deformed blade section; forming the strain displacement relations; and calculating 
stresses and air velocity components in the flap, lag, and torsion, directions. 

These expressions include geometric nonlinearities owing to the assumption of small 
strains and moderate deflections which give rise to numerous higher-order nonlinear 
terms. So an ordering scheme, based on assigning orders of magnitude to the various 
physical parameters is used to reduce the number of terms. The governing equations 
of motion are then obtained using Hamilton's principle. These equations are non- 
linear, partial differential equations in u,v,w, and in torsion deflection. These 
are converted to ordinary differential equations by using Galerkin's method by 
expressing the bending and torsion deflections in terms of generalized coordinates 
and mode shape functions, 


v = Yj rv (*)*j(x) 


J = 1 


w = X 

j=i 


$ = ^ <t> (<|>)6j(x) 
J=1 


(7a) 


and by expressing the induced velocity as 

v. = X + v n + v x cos + v x sin i|* (7b) 

1 U C X 
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where <i» = fit, x = x/R, and are mode shapes, R is the blade radius; and 

Np* an< 3 Nj are the numbers of flap, lag, and torsion modes, respectively, used 

in the analysis. This yields N nonlinear, nonhomogeneous ordinary differential 
equations in terms of modal generalized coordinates V,, W,, and where N is 
the total number of flap, lag, and torsion modes used in the analysis. The equa- 
tions have constant coefficients for mass, damping, and stiffness in the case of 
hovering flight (and periodic coefficients in the case of forward flight). These 
equations are then linearized for small perturbation motions about the equilibrium 
position by expressing the generalized coordinates in terms of equilibrium quanti- 
ties plus small perturbation quantities. 

V j U) = V) + AV j { *> 

WjU) = W 0J + AWj(i|>) (8) 

4>j(<0 = 4> 0 j + A<t>j ( ij») 

Two sets of equations are obtained from this operation: a set of N nonlinear 

algebraic equations in Vqj, Wqj, anc * ^Oj which define the equilibrium position and 
a set of N equations obtained by subtracting the equilibrium equations and dis- 
carding all nonlinear products of perturbation quantities. 

Three more equations are obtained for the dynamic inflow components from rotor 
perturbations in aerodynamic thrust (C T ) and in pitch (C M ) and roll (C L ) moments 
(see dynamic inflow section). The coefficients of these equations are also func- 
tions of equilibrium solution. The quasi-static torsion model equations are 
obtained by dropping torsional inertia and damping terms in the torsion equation. 

The torsion equation is then solved for A<j> . , and the result is substituted into the 
flap-lag equations. This procedure from reference 26 is explained in Appendix B. 

Dynamic Inflow - The total inflow is assumed to consist of a steady trim value 
and a dynamic inflow component as 


Vi = X + V (9) 

At any point on the rotor disk (x,i|>), the dynamic inflow is assumed to have a 
first harmonic representation in terms of inflow distribution (uniform, fore to aft, 
side to side) as 


v = \>q + v q x cos i|> + v g x sin 


( 10 ) 


These velocity components assume the role of degrees of freedom and are 
linearly related to the unsteady thrust, pitching, and rolling moments by first- 
order differential equations: 
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where the subscript aero indicates that C T , C M , and C L are aerodynamic contribu- 
tions only. The nonzero elements of {m} and {£.} are (ref. 17) 


m^ = 8/3" ; m 22 = m^ = -16/45" 

*11 = 1/2v ’ * 22 = a 33 = -2/v 

where v = 2X and 


( 12 ) 


X = ^6 (V 1 + ao n® ^ " ’) 

For rotors with a finite number of blades, however, blade loading rather than 
disk loading must be defined. Therefore, the disk loading is approximated by 

S * ?5 t [J < V* 

k=1 L 0 

C M * Z [ J ( Vk* d * 003 *k ] ( ,3) 

k=1 L*0 J 

C L = wYj [j <C wpV d * sln \ ] 

k=1 L 0 J 

where L is the perturbed aerodynamic force in the flapping direction, 
wp 

The quasi-static inflow formulation, in which the effects of apparent inflow 
mass are neglected, provides a means of approximating the unsteady inflow effects 
without increasing the system dimension. With momentum theory as a basis, quasi- 
static inflow theory leads to the method of an equivalent Lock number y* and drag 
coefficient (refs. 17-20) 


Y 


* 



1 

1 + ao/8v 


] 


(14) 
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This implies that inflow decreases the expected lift and increases the expected 
drag. It should be noted that y* > (C^/a)* are only in the perturbation 
equations and true y, C^^/a are used m calculating the time-independent hovering 
equilibrium displacements. 

Equations through HESL - In this study, a Lagrangian formulation is used to 
obtain the governing nonlinear ordinary differential equations through symbolic 
manipulation. The input to the program are the relations, in alphanumeric format, 
for the position vector, strain expressions, air velocity components, and transfor- 
mation matrices as given by Kaza and Kvaternik (ref. 24). The integration relations 
(if known), differentiation relations, the order of the variables, the ordering 
scheme to be used, and the variables for which coefficients are to be collected, are 
also given as data. The program calculates the strain energy, kinetic energy, and 
the generalized force for the given ordering scheme in generalized coordinates using 
assumed modes as in equation (7a). The order of the variables and the ordering 
scheme as prescribed in reference 15 is used in the present derivation. The pertur- 
bation relations as given in equation (8) are substituted to obtain the steady and 
perturbed terms. 

The program generates both the steady-state and linearized perturbation equa- 
tions and the loading terms necessary for calculating aeroelastic stability and 
response. The thrust, pitch-, and roll-moment equations required in the dynamic 
inflow equations are also generated using the perturbed aerodynamic forces in the 
flap direction, which are generated by the program. Contributions from the equilib- 
rium displacements to the perturbation equations from each mode are included. The 
equations generated are written automatically into FORTRAN subroutines for subse- 
quent numerical calculations. A numerical program subsequently identifies the mass, 
damping, stiffness, and load coefficients for each equation and for each mode. It 
should be noted that for a hingeless rotor, the axial displacement can be solved for 
a priori as a function of flap bending and lag bending. In the present study, 
expressions for axial displacement and axial velocity are taken from reference 15 
and are supplied as data to the program. 

For the results presented here, the program was run on an IBM 360. It took 
approximately 450 sec to derive structural terms and 165 sec for aerodynamic terms 
with inflow (60 sec without inflow) for an elastic flap-lag-torsion analysis. It 
then took 25 sec to write the individual matrix elements into the subroutines for 
subsequent analysis. 

Sample input and output - A brief description of the program input and output 
follows. Table 2 shows the FORTRAN symbol definition used for the original varia- 
bles. Table 3 shows the input to calculate tangential and perpendicular velocities 
using the transformation matrix LAFP(READ MATRIX) and air velocity vector VEL(READ 
MATRIX). By multiplying the two matrices (FOTM MATRIX) with ordering scheme *E2D1, 
the vector AVEL is obtained which gives the components of the velocities in radial, 
tangential, and perpendicular directions. The vector components are redefined as 
expressions by command MATRIX EXPRESSIONS, and since the actual components are a 
negative form of the original expressions, the expressions are negated by calling 
the NEGATE command, giving the actual velocity expressions. 


29 


The following points are to be noted in giving the data to the symbolic pro- 
gram. (1) In hover, the nonlinear rate terms do not contribute to the steady solu- 
tion or the linearized stability solution. Hence, the number of nonlinear terms and 
perturbation relations may be smaller. (2) By redefining the orders of variables at 
proper stages of the input, higher order terms can be retained in lag and torsion 
equations. 

Method of solution - The bending and torsion mode shapes used in equation (7) 
are taken as the uncoupled nonrotating mode shapes of a cantilever beam (ref. 26). 

A total of 15 (N = 15) uncoupled nonrotating modes are used in calculating the 
steady state deflection. The integrals in the problem are evaluated numerically 
with a 16-point Gaussian integration scheme. The steady state equations obtained 
using equation (8) are time-independent in hovering flight. They are nonlinear, 
algebraic equations in Vqj, Wq, , and <t> 0 j • These equations are solved iteratively 
using Brown's algorithm (ref. 27). This algorithm is a modified form of the Newton- 
Raphson method, and the solution is estimated such that the sum of the squares of 
the errors is a minimum. The solution to the linear equations is used as the ini- 
tial estimate in the solution procedure. 

The aeroelastic stability of the blade's motion about the equilibrium position 
is determined by the eigenvalues of the perturbation equations ( A V j , AWj, A<J>j, v 0' 
v c , and v s ). This results in ( (2N+3)x(2N+3) ) matrix for the eigenvalue analysis. 

As explained in reference 15, this can be reduced to a ( (M+3)x(M+3) ) matrix by using 
the lowest M eigenvectors of the free vibration problem. In the present study, 
six (M = 6) eigenvectors, based on the free vibration analysis using 15 (N = 15) 
uncoupled nonrotating modes, are used. The reduced matrices are analogous to stiff- 
ness and damping matrices generated from M-coupled, rotating-blade modes. These 
equations (with M generalized coordinates, say x^, x 2 , ... x M ) are then trans- 
formed to multiblade coordinate equations. The coupled blade and dynamic inflow 
equations are then solved using an eigenanalysis for stability of the rotor. 

Multiblade coordinate transformation - To avoid periodicity, and to have a 
better understanding of the rotor behavior, it is necessary to express the final 
equations in a multiblade coordinate system (a nonrotating frame of reference). As 
pointed out in reference 28, for a constant-coefficient system (as in hover), the 
multiblade coordinate transformation acts only on the degrees of freedom. Let (Xj) 
represent the M-generalized coordinates (obtained above) in the rotating frame, 

{Xj} = {X r X 2 ,X 3 ,X 4 ... X M } (15) 

Then the equations for kth blade including dynamic inflow are 


[M]{X jk } + [C]{x jk } + [k]{X jk } 


+ [T] 


t 


' •- 1 


c — \ 
sin 4>, 


= 0 


(16a) 
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( 1 6b) 



C„, C u , and C. can be expressed in terms of the blade degrees of freedom and 
1 M L 

inflow degrees of freedom as 

b 

G T = 7b Yj ( A 1 j X jk + 6 1 J X jk + C 1j X jk + G 1i u i ) 
k=1 


'M 


-a a 
yb 


I 

k= 1 


( fl 2J x Jk * 8 2J i Jk * c 2J x Jk * G 2lUl )oo3 * k 


(16c) 


-ao 

Yb 


E <» 


3j X jk * S 3J*Jk * C 3J X jk * G 31 u c )si " *k 


k= 1 


where Xj k are the generalized coordinates of kth blade; j runs from 1 to M, the 

number of coupled rotating-blade modes used; i is 1 to 3, the inflow degrees of 
freedom, and U 1 = v Q , U 2 = v Q cos ii» k ; and LU = v g sin . It should be noted that 
matrices [M], [C], [K] , [T], [A], etc., are already reduced to MxM matrices using 
the eigenvectors of a free vibration of coupled rotating beam, as explained in the 
previous section and in reference 15. 

Each generalized coordinate Xj k is written in multiblade coordinates as 

X jk ■ X 0j + X cj 003 *k * X sJ sln *k 


X jk = X 0j * <*cj * X sJ )cos \ * (i sj - X oj )sin *k 


(17) 


x jk = X 0j * <X CJ * 2X SJ - X Cj )C0S *k + <X 3j - 2X Cj - X SJ )S1 " "k 


where Xqj represents collective mode, X c j and X s j represent cyclic modes in the 
nonrotating frame. 
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Let qQ be the vector of the collective modes of all modes; let q c be the 
vector of the longitudinal cyclic modes of all modes; and let q g be the vector of 
the lateral cyclic modes of all modes. Then the multiblade coordinate transforma- 
tion, for a three-bladed rotor is given by 


M 

0 

0 




(18a) 




For a four-bladed rotor, the equation for the differential collective mode is 
not coupled to the equations for the other modes (in hover only). Consequently, the 
same equations apply for four-bladed rotors also, and to rotors with any number of 
blades greater than, or equal to, 3. 

The equations are written for clarity as 

[P] {Y} + [Q] {Y} + [R]{Y} ♦ [T] {u} = 0 (19a) 
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(19b) 


[A]{y} + [B]{Y} + [C]{Y} + [G] {u} = [m]{u} + [a] _1 {u} 

Y = q 0 ,q c ,q s , « = VV v s 

In state vector form, this becomes 

0 0 ° 1 f Y> l 

-P -1 R -P _1 Q -P _1 T m (20) 

m -1 (C - AP -1 R) m _1 (B - AP _1 Q) m _1 (G - t" 1 - AP _1 T)J L'-J 

or 

Z = [A g ]Z 

where [A g ] is a constant-coefficient matrix. The real part of the eigenvalues of 
[ A_ ] gives the modal damping information. 

It can be observed from equations ( 18a) and ( 1 8b ) that the collective mode 
decouples from the progressing and regressing modes. This fact can be used for 
efficient programming in the calculation of eigenvalues. 

Results and Discussion 

To demonstrate the capability of the symbolic processor program discussed, the 
aeroelastic stability of a hingeless rotor system in hover is investigated. As 
schematically shown in figure 1(a), the numerical results presented here are 
obtained directly from the computer-generated state equations. Results are compared 
with previously published damping data to validate the program. Next, the influence 
of the structural model with and without dynamic inflow is studied. In the deriva- 
tion of equations, the orders of magnitude assigned for each parameter is the same 
as that followed in reference 15. The parameters used for numerical study are 

o> v = 0.7, 1.5 , « w = 1.15 , = 3, 5 , K m i /K m2 = 0 ’ 

K /R = 0.025 , c/R = 0.07854 , y = 5.0 , a = 0.05, 0.1 , 
in 

a = 2 ., c do = o.oi, s po = o.o, ayy 2 = i.i5 

Figures 4 and 5 compare the lead-lag damping values obtained in reference 5 
from the equations of Hodges and Dowell (ref. 12), and from the computer-generated 
equations. The results are presented for two torsional frequencies (oj^ = 2.5, 5.0) 
for both soft in-plane (o> v = 0.7, fig. 4) and stiff in-plane (to v = 1.5, fig. 5) 
configurations. The other parameters are u» w = 1.15, R = 1.0, and a - 0.1. It can 
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be seen that each set of equations gives almost identical results for the first mode 
lead-lag damping in the operating range of 0 = 0 to 0.4 rad. The slight differ- 
ences are a result of the following. In deriving the computer-generated equations, 
basic expressions have been taken from Kaza and Kvaternik (ref. 24). The equations 
of references 23 and 24 differ in the structural terms (4>'w'GJ)" in the lag equation 
of motion, -(<t>'v"GJ)' in the flap equation of motion, and (v"w'GJ)' in the torsion 
equation of motion. The aerodynamic ( pac/2) (c/2) [nx( 0 + <|>) + v. ]4> term in the lag 
equation and (pac/2)£2 x v'w' in the flap equation differ between the two aerody- 
namic models. The term 

/>- 

does not appear in the present flap equation which is included in reference 23. The 
kinematics of these two sets of equations are discussed in references 29 and 30. 

The results from these two comparisons provide a validation of the scheme for 
symbolic generation of the equations. Results are now presented for more complex 
rotor modeling problems using the computer-generated equations. 

Elastic blade versus rigid blade (flap-lag model) - The elastic flap-lag-torsion 
equations are reduced to flap-lag equations by dropping all the torsion terms. 

Three modes each for flap and lag are used in the numerical calculations. By care- 
fully studying the rigid blade equations (ref. 13) and the elastic blade flap-lag 
equations (ref. 15), it can be seen that the elastic blade equations have additional 
stiffness contributions in the aerodynamic expressions. This coupling of aerostiff- 
ness is a function of both flap and lag trim values. Figures 6 and 7 show a compar- 
ison of rigid blade and elastic blade results for soft in-plane and stiff in-plane 
rotors, respectively. These figures show the lead-lag mode damping plotted versus 
pitch angle for a rigid blade model, for an elastic blade model, and for an elastic 
blade model with the aerostiffness terms dropped in the elastic blade equations. 

Also shown is lead-lag mode damping plotted for a rigid blade model with no elastic 
coupling effects. It can be observed for both configurations that the elastic blade 
model without the aerostiffness terms shows good correlation with the rigid blade 
equation results, since in the rigid blade model, there is no coupling aerostiffness 
term. The aerostiffness terms seem to have a destabilizing effect on the stabil- 
ity. As noted in reference 12, the rigid blade model with no elastic couplings 
(R = 0.0) yields poor correlation with the elastic blade model for both rotors. It 
should be noted that the damping curves from the present study for the stiff 
in-plane blade in figure 7 are identical to those given in reference 12, again 
validating the symbolic derivation and numerical results of this study. 

The following paragraphs and figures study different blade structural models 
with and without dynamic inflow. The convention used for the figures is one which 
uses solid lines to refer to results without dynamic inflow and dashed lines to show 
results with dynamic inflow. 
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Rigid blade versus elastic blade with dynamic inflow (flap-lag model) - Figure 8 
shows the lead-lag regressing mode damping for the flap-lag model (all torsion terms 
dropped) with dynamic inflow for a rigid blade model, elastic blade model, and an 
elastic blade model with the aerostiffness terms dropped. Again, it can be seen 
that the elastic blade flap-lag equations without aerostiffness terms give good 
agreement with rigid blade results. However, the elastic blade is always less 
stable than is the rigid blade model. It may be noted that the rigid blade results 
with dynamic inflow are identical to those presented in reference 18 for zero elas- 
tic coupling (R = 0.0). 

Elastic flap-lag-torsion with and without dynamic inflow - Figures 9 and 10 show 
the lead-lag regressing mode damping plotted for varying lag frequency for two 
torsional frequencies = 3.0, 5.0). Results for an elastic flap-lag model are 
also shown. By comparing the no dynamic inflow curves with those given in refer- 
ence 15 (figs. 42 and 44), it can be seen that for the results of this study the 
solid-line curves do not cross at the lag frequency that is equal to the matched 
stiffness value (w v = 0.57). This is due to the following reason. In the present 
study, the kinematic pitch-flap and pitch-lag couplings, as mentioned in refer- 
ence 15, are not at zero even for the matched stiffness case. There is a contribu- 
tion to these couplings that is equal to the torsional rigidity times the bending 
trim values owing to the terms mentioned in the second paragraph. 

Figure 9 shows the lead-lag regressing mode damping plotted for zero elastic 
coupling for different torsional frequencies (w^ = 3.0, 5.0), and for a flap-lag 
model, with and without dynamic inflow. It can be observed for the flap-lag model 
that the dynamic inflow effect reduces with increasing lag frequency and is always 
stabilizing. For a torsion frequency of 5.0, the dynamic inflow effects are stabi- 
lizing, and the effect increases with increasing lag frequency. For a torsion 
frequency of 3.0, the dynamic inflow effects are zero or slightly negative, but are 
stabilizing with increasing lag frequency (w v > 0.9). 

Figure 10 shows the lead-lag regressing mode damping plotted for the same 
parameters but for full elastic coupling. It is observed that dynamic inflow again 
increases damping for torsionally rigid blades, and this effect is reduced with 
increased lag frequency as in the previous case. But for torsionally flexible 
blades, dynamic inflow reduced damping for both the torsion frequencies considered, 
except for high values of lag frequency. 

Figure 11 shows the lead-lag regressing mode damping plotted for varying lag 
frequency for three values of the elastic coupling parameter for a blade with tor- 
sion frequency, w. = 5.0. It can be seen that for an elastic coupling value of 0.4, 
the dynamic inflow effects are almost negligible, except for 0.9 < oj v < 1.2. This 
effect of the elastic coupling parameter R on damping with dynamic inflow is 
further investigated in figures 12 and 13. Here the lag regressing and flap 
regressing mode damping are plotted with varying elastic coupling for a flap-lag 
model and for a flap-lag-torsion model with torsional frequencies of 3.0 and 5.0. 

The lead-lag frequency is w v = 1.4. Figure 12 shows the lead-lag regressing mode 
damping. It can be seen that for torsionally flexible blades, depending on the 
value of the coupling parameter R, dynamic inflow can be stabilizing (R < 0.3) or 
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destabilizing (R > 0.4). Consequently, there exists an elastic coupling value for 
which dynamic inflow effects on lead-lag damping are negligible. This can also be 
interpreted as showing that for torsionally soft blades, the dynamic inflow reduces 
the magnitude of the lag mode damping regardless of sign. The value of R for 
which dynamic inflow effects are negligible is where the lag damping is 
approximately zero. Noting that the points at which lag damping with and without 
dynamic inflow are zero and almost identical , it may be stated that dynamic inflow 
does not have a significant impact on the stability boundary. These results are 
different from those for the elastic flap-lag blade for which dynamic inflow effects 
are always stabilizing for the lead-lag mode. In figure 13, the flap regressing 
mode is reduced with dynamic inflow for all values of elastic coupling parameter 
(0.0 < R < 1.0), as expected. 

The influence of the dynamic inflow model is further investigated in 
figure 14. The lead-lag regressing mode damping is plotted for varying elastic 
coupling values for a dynamic inflow model and for a quasi-static inflow model. 

Also shown are the damping values obtained by taking 9 = \ (eq. (12)), as in refer- 
ence 21. The dynamic inflow models of references 17 and 21 give similar results, 
particularly at high elastic coupling parameters. The quasi-static model is a 
reasonable approximation for this case. 

Figure 15 shows the lead-lag regressing mode damping plotted from a quasi- 
static torsion model with and without dynamic inflow (w^ = 5.0). It also shows a 
comparison with a coupled flap-lag-torsion analysis. The quasi-static torsion model 
has been obtained by dropping torsional inertia and by damping the terms in the 
torsion equation. The torsion equation is solved for the perturbed torsion deflec- 
tion (A<t> <) , and this is substituted in flap-lag equations. It can be seen that the 
quasi-static torsion model exhibits the same qualitative behavior as the flap-lag- 
torsion model does, although the in-plane damping is uniformly reduced for all 
values of R. It may also be observed that dynamic inflow again reduces damping for 
elastic parameter values greater than 0.4. 

As was explained in reference 15, the quasi-static torsion model attempts to 
represent the torsion effect as equivalent kinematic pitch-flap and pitch-lag cou- 
plings. These couplings depend on the blade's torsion frequency, the elastic cou- 
pling parameter, equilibrium bending deflections, and the difference in lead-lag and 
flap bending stiffness. In the present study it has been shown earlier that these 
couplings also depend on torsional rigidity. For clarity, a rigid blade model with 
pitch-flap and pitch-lag couplings is studied. Results are presented with and 
without dynamic inflow for a varying elastic coupling parameter. Figure 16 shows 
the lead-lag regressing mode damping versus the elastic coupling parameter from a 
rigid blade flap-lag model having pitch-flap and pitch-lag couplings with and with- 
out dynamic inflow. Results from the quasi-static torsion model are also plotted 
for comparison. For a given flap, lag, and torsion frequency, the values of pitch- 
flap and pitch-lag couplings are calculated from equation (40) of reference 15. For 
comparison, the trim-bending values obtained from the quasi-static torsion model are 
used to calculate these couplings. For u> v = 1.4, u> w = 1.15, = 5.0, and 

9 = 0.3, it is observed that the value of pitch-flap coupling varies from +0.125 
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for R = 0.0 to -0.414 for R = 1.0. The pitch-lag coupling varies from -1.4602 
for R = 0.0 to -1.497 for R = 1.0. Comparing the curves for both the models with 
no dynamic inflow, it can be seen that the rigid blade model yields conservative 
results. With dynamic inflow, both models qualitatively exhibit the same trend, 
namely, that a reduction in damping results with increasing elastic coupling for 
R > 0.5. It may be concluded from figures 15 and 16 that a rigid blade flap-lag 
model with pitch- flap and pitch-lag couplings, rather than a model with coupled 
flap-lag-quasi-static torsion for full flap-lag-torsion, may in some cases be suffi- 
cient to yield qualitative trends. 

Summary of Hover Study 

The combined effect of blade torsion and dynamic inflow on the aeroelastic 
stability of an elastic rotor blade in hover has been studied. The governing equa- 
tions of motion of the elastic blade with flap-lag-torsion degrees of freedom and 
the dynamic inflow equations are derived using a symbolic processor written in 
FORTRAN. The blade and the dynamic inflow equations are converted to equations in a 
multiblade coordinate system by rearranging the coefficients of the equations. 
Conclusions drawn from this study are presented in the last section of this report. 

Forward Flight 

In forward flight the rotor does not have an axisymmetry flow environment as it 
does in hover. The velocities of the rotor blade not only depend on the blade 
motion and on rotational speed but also on the on-coming air speed. These addi- 
tional components of air velocity are shown in figure 17. The advancing blade has a 
velocity relative to the air higher than the rotational velocity, whereas the 
retreating blade has a lower velocity relative to the air. The aerodynamic forces 
vary periodically as the blade rotates. This causes asymmetry of the loads. The 
rotor blade loading and motion are periodic with a fundamental frequency equal to 
the rotor speed fl. The analysis is more complicated because of the dependence of 
the loads and motion on the azimuth angle U> . A phenomenon which can be introduced 
by forward flight is the reverse-flow region, an area on the retreating side of the 
rotor disk where the velocity relative to the blade is directed from the trailing 
edge to the leading edge. In general, the reverse-flow region is defined as the 
area of the disk where the tangential velocity is less than zero, which has the 
boundary x + p sin = 0. This yields a circular region with a diameter p, cen- 
tered at x = p/2, and ip = 270 deg on the retreating side. Aeroelastic studies in 
forward flight are consequently more complex than are hover studies. The governing 
equations are lengthy and have periodic coefficients. However, as in hover, 
reliable solutions can be obtained by linearizing the nonlinear equations about a 
steady-state equilibrium position. In forward flight this position is a time- 
dependent periodic solution. Calculation of this time-dependent equilibrium posi- 
tion is inherently coupled with the trim state of the complete helicopter in forward 
flight. 

Formulation of forward flight stability analysis - The blade model remains the 
same as described for hover and shown in figure 2. The blades undergo axial (u), 
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lagwise benging (v), and flapwise bending (w) and torsion deflection ( <t» ) - The 
aerodynamic forces will have additional contributions from the airspeed. These 
contributions have been shown in figure 17. The general formulation of the problem 
follows the same procedure as the one for hover. However, in forward flight the 
implementation of Galerkin's method requires the rotating modes to be calculated at 
every instant of time for a given pitch. It is proposed here to derive the govern- 
ing equations of motion using the uncoupled rotating modes evaluated at zero 
pitch. The elastic degrees of freedom v, w, and <j> are expressed in terms of the 
generalized coordinates and mode shape functions as 


w = 2 ^ Rw^ilOn^x) 
i=1 


v = £ Rv^ipk^x) (23) 

i = 1 


<> = 2^ 

i= 1 

where ip = fit, x = x/R, r^, ? it and 9. are rotating mode shapes; N p , N L , and N T 
are the number of flap, lag, and torsion modes used in the analysis. These mode 
shapes are obtained from five nonrotating modes for each degree of freedom. In the 
actual derivation of the equations through HESL, the integral quantities are evalu- 
ated by the numerical program using these modes. Two rotating modes are used in 
this study for each flap, lag, and torsion degree of freedom. 

Equations through HESL - The derivation of the governing equations of motion for 
a single blade using HESL is the same as described for hovering flight. The only 
change in data input is the air velocity components which are caused by the air 
speed. The data for the position vector, strain expressions, air velocity compo- 
nents, and transformation matrices from reference 24 are given as the basic data. 

The displacements as given in terms of the generalized coordinates (eq. (23)) are 
substituted and the equations are derived using a Lagrangian approach. The per- 
turbed equations are obtained by substituting the relations given by equation (8). 
The nonlinear single blade equations are written to a separate file, to be used in 
the time-dependent equilibrium position calculation. The perturbed equations are 
similarly written to a file to be used in a multiblade coordinate transformation. 

The relations required in dynamic inflow calculations, thrust (C^), pitching moment 
(C M ), and rolling moment (C L ) are derived using the perturbed aerodynamic force in 
the flap direction. The equations are lengthy because of the contributions from 
airspeed. The program generates the steady and perturbed equations in a single 
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operation and outputs them individually. This is convenient in the case of forward 
flight not only because of the large number of terms present in each equation, but 
also for the different analysis processes that are required for the steady and 
perturbed sets of equations. 

In the present study, the order of the variables and the ordering scheme used 
are same as those followed in reference 15. All the o(e^) terms, when compared 
to o(1), except those that contribute to ^ead-lag and torsion damping are 
neglected. Nonlinear rate products (vw, v , etc.) are retained since they contrib- 
ute to the linearized stability analysis. Although any general ordering scheme 
could have been used to obtain the final equations of motion, this ordering scheme 
is considered representative and adequate for demonstrating the capability of the 
symbolic analysis process. For the results presented here, it took approximately 
300 sec to symbolically derive both the structural and aerodynamic equations on a 
VAX 11/780 computer. It should be noted that for a hingeless rotor, the axial 
displacement can be solved for a priori as a function of flap and lag bending. In 
this report, expressions for axial displacement and axial velocity are taken from 
reference 15 and supplied as data to the program. 

Dynamic inflow modeling - The dynamic inflow equations are related to the blade 
degrees of freedom through the variations in thrust, pitching and rolling moments 



The elements of (m) and (a) define the various dynamic inflow models that can be 
included in an analysis. Reference 19 presents a hierarchy of models, having dif- 
ferent elements for (m) and (a) from actuator disk theory in forward flight. The 
elements of (a) depend on the wake skew angle at the rotor 
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. -IX 
tan — 
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(25) 


where X is the steady inflow. Of the 13 models presented in reference 19, the 
partially constrained model gave good results. In the present report, this par- 
tially constrained theory is used to obtain the dynamic inflow results. The ele- 
ments of (m) and (a) are given by 
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Trim and periodic equilibrium solution - The nonlinear periodic coefficient 
equations obtained earlier can be solved for the periodic response in the time 
domain using a Floquet method or in the frequency domain using a harmonic balance 
method. Either will yield the time-dependent equilibrium position about which the 
nonlinear equations can be linearized for an eigensolution. In forward flight, this 
equilibrium position is coupled with the entire trim state of the helicopter. The 
trim state is the steady state condition achieved by the system as time increases 
without bound with the controls fixed and no external output present. Calculation 
of a trim position requires establishing the control settings for a given flight 
condition. The control settings are collective pitch, longitudinal and lateral 
cyclic pitch, and the rotor shaft angle of attack. The induced velocity which 
depends on the generated thrust and advance ratio is also calculated. In this 
report, the harmonic analysis method coupled with an iteration on the trim state is 
employed to calculate the equilibrium position and the trim settings. This numeri- 
cal procedure, taken from reference 26 is as follows. An iterative inner loop in 
which the solution for the periodic motion is obtained with fixed values of the trim 
variables is nested within an outer loop in which the solution for the trim varia- 
bles is obtained. The rotorcraft motion is solved for the periodic motion by the 
harmonic analysis method which directly calculates the harmonics of a Fourier series 
representation of the motion. The procedure advances the rotor around the azimuth 
calculating the forcing functions in the time domain and then updating the harmonics 
of the response. The forces and moments acting on the rotor are calculated from 
this response, and the controls are adjusted until the equilibrium of forces and 
moments that are required for the specified operating state is achieved. 

For the harmonic analysis method, the governing equations of motion are written 
with all the time-dependent and nonlinear terms as a forcing function, as 

[M] {x} + [C] {X} + [K] {x} = F(X,X, 1 ») (27) 

where matrices M, C, and K are the constant parts of the mass, damping and stiff- 
ness matrices, and where X is the vector of degrees of freedom. The function F 
is evaluated at J points around the rotor azimuth 

F = F(* ) (28) 
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and the harmonics of a complex Fourier series representation of F are given by 


F 


n 



-inip. 
e J 


(29) 


Then the nth harmonic of the motion is given by 

X = H -1 F (30) 

n n n 

p 

where H = K - n M + iCn. The iterative solution proceeds as follows. At a given 
azimuth ipj, the blade motion is calculated using the current estimates of the 
harmonics 
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where L is the number of terms used in the expansion of X. The forcing function 
is evaluated with this motion. Then the change in the harmonics caused by the 
change in the forcing function is calculated and added to the harmonics calculated 
in the last revolution. After every revolution the old and new harmonics are 
checked until convergence is achieved. 

After obtaining the harmonics of the blade motion, the equilibrium of the rotor 
forces and moments is checked. If overall equilibrium is not satisfied, the trim 
settings are changed and the procedure is repeated until equilibrium is obtained. 

The following assumptions are made in arriving at the equilibrium of the forces and 
moments. The helicopter is in straight and level steady flight, the rotor hub 
moments are trimmed to zero with fuselage and other aerodynamic contributions 
neglected. Then the equilibrium of forces and moments are given by 

C,p cos a + sin a = 

-C T sin a + cos a = -C 

C = C =0.0 
mx my 

_ 2 

where = (1/2)fy and a is the shaft angle of attack. 


1 


( 32 ) 
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In calculating thrust, horizontal force, and the hub moments, the steady inflow 
appears as a parameter which in turn depends on the rotor thrust and shaft angle of 
attack. In this study, the study inflow is taken as an equation of constraint and 
solved along with the four equations above. 
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y tan a 
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(33) 


The increments in the trim settings are calculated using a modified Newton Raphson 
method (ref. 30). If v is the control variable and M is the target to be 
achieved, then a first-order approximation of M(v) is 


"target = "n.1 = "n * <Vl ' v n> 


or v n+| = v + D~ (^target ' "n )f <31,) 


where v and v -| indicate the current and new estimates of v respectively and 
f < 1 is included to avoid overshoot oscillations in the trim iteration by reducing 
the step size. The partial derivative matrix D is 
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(35) 


where v- is the ith control variable and Av^ is its increment. The conver- 
gence is checked when the tolerence level as specified by the parameter e is 
met. For more details on this procedure see reference 30. In the present study, 
all the degrees of freedom, blade flap, lead-lag and torsion degrees of freedom are 
used in calculating the trim state and periodic response. 

The periodic response can also be obtained using Floquet methods (ref. 16). A 
review of these methods to obtain the initial conditions, forced response and sta- 
bility data is given in reference 31. Explicit harmonic balance equations can also 
be derived symbolically using the program HESL. This requires the definition of the 
harmonics for each degree of freedom and a table of trigonometric relations giving 
the products of sines and cosines as sums of sines and cosines. Then the SUBSTITUTE 
TABLE INTO RELATIONS command substitutes these relations into the governing equa- 
tions. A numerical program collects the coefficients of sines and cosines and forms 
the harmonic balance equations. However, this procedure requires modification of 
the data for each degree of freedom when a larger number of harmonics is involved. 
Hence, HESL is not convenient to use to derive explicit harmonic balance equations 
and numerical method as explained above is employed in this study. 

Aeroelastic stability solution (single blade) - Once the time dependent equilib- 
rium position is determined, the nonlinear equations are perturbed about this equi- 
librium position, equation (8). As mentioned earlier, the symbolic program gener- 
ates the perturbed equations and writes them into subroutines. It should be noted 
that squares of the perturbation quantities are typically neglected by the symbolic 
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program when employing a realistic ordering scheme. The final system of equations 
for stability is 

[P] {Aq} + [Q] {Aq} + [R]{Aq} = 0 (36) 

or 

Z =[AjZ 

The stability of this linearized system is determined from Floquet theory by 
evaluating the characteristic exponents of [A g ] 

X k = 5 k + iu, k (37) 

The mass, damping and stiffness terms are identified by a numerical program for each 
time step. The linearized system is stable when < 0. 

Multiblade coordinate transformation - It is necessary to convert the equations 
into a fixed coordinate system to provide a better understanding and to include 
dynamic inflow effects which are referenced to the fixed system. The multiblade 
(Fourier) coordinate transformation is a linear transformation of the degrees of 
freedom from the rotating to nonrotating frame. In the case of hovering flight, 
this transformation has been easily performed by rearranging the coefficients of the 
equations, since the coefficients are constant. For forward flight the coefficients 
are periodic. Hence the multiblade coordinate transformation is more complicated. 
Let X be a generalized degree of freedom (dimensionless) in the rotating frame for 
the ith blade. Then for a three-bladed rotor, the relations 

X. = X rt + X cos i|>. + X sin 4j. 
i 0 c is l 

Xj = X 0 ♦ <X C ♦ X s )cos * i . (X s - X c )sln (ij (38) 

h -- * 0 * (X c * (X c * 2i s ' X c )cos *i + <*s ' 2X c - V sln *1 

give the ith generalized blade degree of freedom using multiblade coordinates in 
the nonrotating frame. The variables X Q , X c , and X g are the rotor degrees of 
freedom, and they describe the motion of the rotor as a whole in the nonrotating 
frame. 

The Multiblade Coordinate Transformation involves the following steps 
(ref. 29): (1) expansion of each degree of freedom into multiblade coordinates; 

(2) multiplication of the resulting expression with multiblade functions 1, cos 4> , 
sin 4>, cos 24 ), sin 24 ), etc., depending on the number of blades; (3) replacement of 
products of sines and cosines as sums of sines and cosines using trigonometric 
relations; and (4) deletion of terms which are not periodic multiples of the number 
of blades (summation rules). Using the symbolic program, this is achieved as fol- 
lows. The multiblade expansion of each degree of freedom and its time derivative is 
given as a table of relations (eq. (38)). The multiblade functions 1, cos 4>» 
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sin i|>, cos 2i|», sin 2 ip , etc., required in transforming the equations, are read as 
data (or can be generated within the program). The trigonometric relations giving 
the product of sines and cosines as sums of sines and cosines are given as a table 
of relations. This requires that the equation derivation consist of several runs 
until all the required relations are included in the table of relations without 
exception. Then the command PERFORM MULTIBLADE TRANSFORMATION multiplies the equa- 
tion with each of the multiblade functions, substitutes the multiblade expansion for 
each degree of freedom, substitutes the trigonometric relations (from the given 
tables of relations), and checks for the multiples of the blade harmonics. Only 
terms which contain a periodic multiple of number of blades are retained. The 
interface routines separate the terms into constant and periodic parts (coefficient 
of each harmonic separately), and writes them into subroutines. In this manner, the 
constant coefficient approximation is easily done. In this study the perturbed 
governing equations of motion and perturbations in thrust and in moment equations 
are converted into multiblade coordinate equations. The multiblade solution was 
subsequently checked for accuracy by comparing results with a single blade solution 
without dynamic inflow. It should be noted that the trim value harmonics entering 
as nonlinear contributions should be defined as symbolic data. The input data 
increase as more nonlinear terms are taken into account. These data are in addition 
to the data given for the multiblade expansion of the degrees of freedom. However, 
the output may be smaller since only terms which are multiples of the number of 
blades are retained. Since this is a feasibility study to obtain explicit multi- 
blade equations using a symbolic program in FORTRAN, the nonlinear quantities are 
assumed to provide only first harmonic forcing contributions. For the results 
presented here, the program was run on a VAX 11/780 computer. For each equation, it 
took approximately 250 sec to derive the multiblade equations and 120 sec to write 
these into subroutines for numerical analysis. 

From the amount of input required to perform a symbolic formulation of the 
harmonic balance and multiblade equations, the program HESL is convenient to explic- 
itly consider the symbolically derived equations only if the number of degrees of 
freedom are small. As pointed out in reference 29, numerical schemes are better 
suited for general models so that the harmonic balance equations and multiblade 
equations can be efficiently obtained after obtaining the steady and perturbed 
equations from the symbolic program. 

The symbolic program separates the terms containing the periodic variable 
cos Nt and sin Nt and writes the equation as 

A(t) = A_. + A M cos Nt + B m sin Nt (39) 

UN N 

In subroutine form they are referred to as A(1), A(2), and A(3). This allows for 
direct elimination of the matrices A(2) and A(3) for a simple constant coefficient 
analysis. 

Aeroelastic stability solution (multiblade and dynamic inflow )- The final gov- 
erning equations of motion can be written as 
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[P ] {q} + [Q] {q } + [ R ] {q} + [T]{u} = 0 


(40a) 


for the blade equations, and 

[A] {q} + [B] {q} + [C]{q} + [G]{u} = m{u} + U] -1 {u} (40b) 


for the dynamic inflow equations, where {q} is {qQ, q c , q g } , qQ is vector of all 

collective modes, q c is vector of all lateral cyclic modes, q g is vector of all 

longitudinal cyclic modes, and {u} is (v n ,v ,v ). 

u c s 


Equations (40a) and (40b) can be combined as 


where 
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(41) 


The final stability equations in state vector form are 

{Y} = [A S ]{Y} (42) 

where 


{*} 



The stability results are obtained by calculating the eigenvalues of [A s ], in a 
similar manner to the single blade case. 


The size of the state matrix depends on the number of modes and blades. For 
the flap-lag-torsion model with two modes each, the size of the state matrix is 
36 x 36 without dynamic inflow and 39 * 39 with dynamic inflow. The corresponding 
values for a flap-lag model are 24 x 24 and 27 x 27, respectively. 


Results and discussion - Results are presented for a uniform blade with zero 
built-in twist, zero blade offsets, and zero precone angle. Reversed flow effects 
are neglected. A three-bladed rotor is considered. Two rotating modes for each 
flap, lag, and torsion degrees of freedom are used in the calculations. These modes 
are calculated at zero pitch and are obtained from five nonrotating modes. Results 
for both a single blade solution and a multiblade solution are presented for differ- 
ent blade structural models. All results are for a propulsive trim condition, 
calculated for a specified weight coefficient of C w = 0.007 and an equivalent drag 
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area D/q = f = 0.012(irR ), where D is the drag force and q is the dynamic 
pressure. 

In the derivation of the equations, the orders of magnitude assigned for each 
parameter are the same as that which followed in reference 15. The other parameters 
used for the numerical study are 

u> v = 0.7, 1.4 , u w = 1.15 , = 3.0 , c/R = 0.07854 , Y = 5.0 , 

o = 0.1 , a = 2 tt , C , = 0.01 , 8 = 0.0 , K ,/K _ = 0.0 , 

’ ’ d ’pc ’ ml m2 ’ 

o r 

V R = 0.025, (VK/ = 1.5 

Lead-lag damping values (real part of the characteristic exponent) are pre- 
sented for a soft in-plane and stiff in-plane rotor with and without dynamic 
inflow. The results present an investigation of: (1) the effect of degrees of 
freedom used in the trim analysis on the lead-lag damping, (2) the effect of using 
only one torsion mode, (3) the inclusion of a dynamic inflow model, and (4) the 
difference between periodic and a constant coefficient approximation (where all 
time-dependent coefficients are neglected). 

Influence of trim analysis. The effect of the number of degrees of freedom 
used in the trim analysis on the lead-lag damping is investigated in figures 18 
and 19. Figure 18 shows the lead-lag damping plotted versus advance ratio for a 
soft in-plane rotor (w v = 0.7). It can be seen that a flap-lag-torsion stability 
analysis from a flap trim analysis underpredicts the lead-lag damping. The second 
mode shows the same trend with the difference in predicted damping increasing with 
advance ratio. Figure 19 shows the lead-lag damping plotted for a stiff in-plane 
rotor (w v = 1.4) as a function of advance ratio. The results also show an increase 
in damping when a flap-lag-torsion trim analysis is used. It is also noted that at 
an advance ratio of 0.37 < y < 0.41, the roots separate, and one root becomes less 
stable whereas the other becomes more stable. The damping does reduce as the 
advance ratio is increased beyond y = 0.37. The second mode remains stable at all 
advance ratios considered. 

The increase in damping observed above for both soft in-plane and stiff 
in-plane rotors is due to the different time-dependent equilibrium positions that 
are used to linearize the full nonlinear equations. A full flap-lag-torsion trim 
analysis is consistent in that the blade model has the same degree of complexity in 
both the trim and in the stability analysis. It should be noted qualitatively that 
the same type of trend is observed in reference 1 6 , for both soft in-plane and stiff 
in-plane rotors. This verifies the symbolic and numerical programs for the single 
blade results and forms the basis to check the symbolically derived explicit multi- 
blade equations (and numerical results) presented in the next section. 

Influence of elastic torsion model. Figures 20 and 21 present the lead-lag 
damping plotted versus advance ratio from a flap-lag model, flap-lag-torsion model 
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(two inodes each), and flap-lag-torsion model with only one torsion mode, for a soft 
in-plane rotor (w v = 0.7) using a single blade analysis. Figure 20 presents the 
damping results for full elastic coupling (R = 1.0). It can be seen that the flap- 
lag model underpredicts the lead-lag damping. The model with only one torsion mode 
increases the damping above the model with two modes each. Figure 21 shows the 
lead-lag damping value plotted for zero elastic coupling (R = 0.0). The damping 
levels are much reduced when they are compared with the full elastic coupling 
case. However, the flap-lag model is again lowest damped. 

Lead-lag damping is plotted for a stiff in-plane rotor (co y = 1.4) with varying 
advance ratio in figures 22 and 23. Figure 22 presents the damping results for full 
elastic coupling. The same trend that was observed for the case of a soft in-plane 
rotor (fig. 6) exists. Note that root splitting for high-advance ratios occurs even 
when only one torsion mode is used. 

Figure 23 presents the lead-lag damping for increasing advance ratio for a 
stiff in-plane rotor for a zero elastic coupling parameter and as a result, the 
rotor is unstable. However, a flap-lag model predicts a stable system. This demon- 
strates the importance of elastic blade torsion in the analysis. 

Lead-lag damping is plotted in figure 24 for an advance ratio of y = 0.25 
while varying the elastic coupling parameter for a stiff in-plane rotor. Here a 
flap-lag model predicts positive damping for all values of R, whereas for a flap- 
lag-torsion model, the damping varies with elastic coupling parameter, increasing 
with elastic coupling parameter. 

Multiblade equation results. The following figures present the lead-lag 
regressing mode damping results obtained using multiblade equations. The multiblade 
equations are explicitly derived using the symbolic program. This required explicit 
definition of all nonlinear contributions and degrees of freedom in terms of their 
harmonics. This significantly increased the amount of data that is required by the 
symbolic program. Since this is a feasibility study on the use of symbolic programs 
in FORTRAN, only first harmonics were considered in the nonlinear contributions. 
Consequently, damping data determined from the multiblade equations may slightly 
differ from the single-blade solution. Additionally, the multiblade results are 
obtained by retaining only one torsion mode; however, the nonlinear contribution 
from both torsion modes is used. This reduces the time step and time for integra- 
tion in performing the Floquet-analysis . 

The damping values are first checked with those obtained from a single-blade 
solution obtained earlier to validate the multiblade equation derivation process. 

It was found that the approximation cos 9 = 1.0, used in deriving the explicit 
multiblade equations, will predict slightly higher (yet less than 2 %) damping for 
stiff in-plane rotors with a full elastic coupling parameter, since this approxima- 
tion most greatly affects the coupling elements. For all other values of the elas- 
tic coupling parameter, this approximation does not affect the damping value. Where 
required for comparison, the single-blade damping values are recalculated using this 
approximation, to avoid the rederivation of the multiblade equations. 
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Figures 25 and 26 show the lead-lag regressing mode damping plotted for varying 
advance ratio with and without dynamic inflow from a flap-lag-torsion and flap-lag 
model for a soft in-plane rotor. Figure 25 shows the damping for full elastic 
coupling (R = 1.0). For the flap-lag-torsion model, the dynamic inflow reduces the 
damping at all advance ratios. Its effect is almost negligible at advance ratios 
0.15 to 0.25. For the flap-lag model, the dynamic inflow increases the damping up 
to an advance ratio of 0.33, it then reduces the damping. Figure 26 presents the 
damping results for the zero elastic coupling (R = 0.0). It is seen that for the 
flap-lag-torsion model, the dynamic inflow again reduces the damping in hover. Yet 
at intermediate advance ratios, dynamic inflow increases the damping, and at higher 
advance ratios, it once again reduces the damping. For the flap-lag model, the 
dynamic inflow increases lead-lag damping for all the advance ratios that are con- 
sidered. This is consistent with previous studies, for example, as in reference 18. 

The lead-lag regressing mode damping is plotted for a stiff in-plane rotor for 

varying advance ratio in figures 27 and 28. Figure 27 is for a rotor with full 
elastic coupling. For a flap-lag-torsion model, the dynamic inflow reduces the 
damping up to an advance ratio of 0.41. For y > 0.41, this model shows a slightly 
increased damping value. The flap-lag model with dynamic inflow shows a small 
decrease in damping, yet this damping increment gets smaller with higher advance 
ratios. Figure 28 is for a rotor with zero elastic coupling. For this configura- 
tion, the dynamic inflow increases damping for all advance ratios. Consequently, 
both flap-lag-torsion and flap-lag models show the same trend. 

Figure 29 shows the lead-lag regressing mode damping as it is plotted for a 
stiff in-plane rotor at an advance ratio of 0.25 while varying the elastic coupling 
parameter. For the flap-lag-torsion model, dynamic inflow reduces the damping for 
R > 0.3; however, it increases the available damping for R < 0.3. This increase is 

not sufficient to stabilize the in-plane mode. With the flap-lag model, dynamic 

inflow shows an increase in damping for all values of elastic coupling. 

Constant coefficient approximation. The effect of a constant coefficient 
approximation is studied in figure 30, where the real part of the exponent is 
plotted for a stiff in-plane rotor with full elastic coupling which showed a split- 
ting of the roots when analyzed with a full periodic coefficient. The constant 
coefficient approximation does not show this splitting, since the frequencies are 
away from the real axis. For this analysis, the regressing and collective modes do 
predict the same damping trend with advance ratio as shown by a full Floquet analy- 
sis. However, the progressing mode shows poor agreement between a constant 
coefficient approximation analysis and a Floquet theory analysis. This is because 
the constant coefficient approximation can only be good for low frequency modes. 


III. CONCLUSIONS 


A symbolic manipulation program written in FORTRAN is used to derive the aero- 
elastic analysis equations of an elastic blade with flap-lag-torsion degrees of 
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freedom in hover and in forward flight. A study is made of the feasibility of using 
the program to obtain explicit equations in a harmonic balance method and to obtain 
multiblade equations. Numerical results are presented for a flap-lag-torsion model 
and for a flap-lag model with and without dynamic inflow for a propulsive-trimmed 
rotor. Both soft in-plane and stiff in-plane rotors are considered. 

The following conclusions are drawn from using a symbolic program: 

1. The program HESL, written in FORTRAN, can be conveniently used to derive 
elastic blade equations. 

2. The process of aeroelastic stability analysis starting from equation deri- 
vation to numerical calculation can be automated. 

3. With the present program capability, the amount of data to the symbolic 
program increases rapidly with degrees of freedom and the number of harmonics that 
have been analyzed. 

4. In deriving the explicit harmonic balance equations and multiblade equa- 
tions, the following should be noted: 

(a) To obtain the harmonic balance equations, a numerical method is sug- 
gested since an arbitrary number of harmonics can be used without increasing the 
input data to the symbolic program. 

(b) To obtain the multiblade equations, the perturbed equations in their 
Fourier series form are derived using the symbolic program. Then the multiblade 
equations themselves are obtained numerically. 

It is recommended that a selective combination of the symbolic and numerical 
programs is required for the derivation and numerical study to be efficiently con- 
ducted. 

The following conclusions are drawn from the numerical study performed for a 
hovering rotor. 

1. The unsteady aerodynamic effects on lag regressing mode damping, modeled as 
dynamic inflow, are dependent on the elastic coupling parameter for torsionally 
flexible blades. For zero elastic coupling, dynamic inflow increased the lag 
regressing mode damping. For other elastic coupling parameter values, dynamic 
inflow either increased or decreased lag regressing mode damping. 

2. Use of a quasi-static torsion blade model is conservative, and yields the 
same trend with varying elastic coupling parameter as flap-lag-torsion model. 

3. When compared with a flap-lag model, it is seen that the relationship 
between dynamic inflow and lag regressing mode damping can be a result of the effect 
of pitch-flap and pitch-lag couplings, which depend on the elastic coupling 
parameter. 
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The following conclusions are drawn from the numerical study of a single blade 
solution in forward flight. 

1. A flap-lag-torsion stability analysis from a trim procedure where only the 
flap degree of freedom is used underpredicts the lead-lag damping. 

2. In the case of stiff in-plane rotors, high-forward flight speed is desta- 
bilizing. At high advance ratios, a splitting of the roots is encountered which 
yields two real values at the same frequency. 

3. Using only one torsion mode increases the damping value from the flap-lag 
structural model. 

4. The damping values for a stiff in-plane rotor are sensitive to elastic 
coupling parameter. 

The following conclusions are drawn from the numerical study of a multiblade 
solution with dynamic inflow in forward flight. 

1. For a flap-lag model, for both soft in-plane and stiff in-plane rotors, 
with zero elastic coupling, the dynamic inflow increased damping at all advance 
ratios considered; with full elastic coupling the dynamic inflow increased the 
damping at low-advance ratios, but reduced damping at high-advance ratios. 

2. For a flap-lag-torsion model, dynamic inflow slightly reduced lead-lag 
regressing mode damping for all advance ratios. The same trend is observed for both 
soft in-plane and stiff in-plane rotors. 

3- For a given advance ratio, the variation of damping with an elastic cou- 
pling parameter for a stiff in-plane rotor showed the same trend as it did for the 
hover case. 

4. The constant coefficient approximation for the stiff in-plane rotor does 
not show the splitting of the rotors, since the frequency of the lag mode is away 
from the real axis. 
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APPENDIX A 


This appendix shows the input data (Table A1) and the output data (Table A2) 
from the symbolic program to derive flap equations for a rigid rotor blade. The 
problem definition starts with a flap-lag transformation but is subsequently reduced 
to flap model. This simplified model is used to clarify the program aspects rather 
than the modeling aspects. Even though the data is given at the matrix level to 
show that the symbolic program is capable of handling matrices, in this example it 
is simple to use the definition of the position vector directly as input. Both the 
single blade and multiblade equations for a three-bladed rotor are obtained from the 
symbolic program. 

Table A3 gives the FORTRAN symbols used to define the original variables. 


51 


TABLE A1 .- INPUT TO THE PROGRAM HESL TO OBTAIN LINEARIZED FLAP EQUATION AND 
MULTIBLADE EQUATIONS FOR THREE-BLADED ROTORS 


*§# 

READ MATRIX 
ITFLP0303 
01 

1.0 COSB 

01 
0.0 
01 

-1.0 SINB 

01 
0.0 
01 
1.0 
01 
0.0 
01 

1.0 SINB 

01 
0.0 
01 

1.0 COSB 

READ MATRIX 
ITLAG0303 
01 


1.0 

01 

- 1.0 
01 
0.0 
01 
1.0 
01 
1.0 
01 
0.0 
01 
0.0 
01 
0.0 
01 
1.0 

READ MATRIX 
#P0VE0301 
01 
1.0 
01 
0.0 
01 


COSZ 


SINZ 


SINZ 


COSZ 


RB 
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TABLE A1 . - CONTINUED 


READ DIFERENTIATION TABLE 

0DERV1O 


BT BTOI 

1.0 


BTD BTD01 

1.0 


BT C0SB01 

-1.0 

SINB 

BT SINB01 

1 .0 

COSB 

TAU BTOI 

1.0 

BTD 

TAU BTD01 

1 .0 

BTDD 

TAU C0SB01 

-1.0 

SINB BTD 

TAU SINB01 

1.0 

COSB BTD 

TAU CSCY01 

-1.0 

SNCY 

TAU SNCY01 

1.0 

CSCY 

READ TABLE 

FOR SUBSTITUTION 

&SUPR02 


01 SINZ01 


0.0 


01 COSZOI 


1.0 


READ TABLE 

FOR SUBSTITUTION 

0DUMYO1O1 


02 SINB02 


1.0 


-1.0 

COSB COSB 

READ TABLE 

FOR SUBSTITUTION 

0DUM1O1 


01 ABCD01 


1.0 


READ TABLE 

FOR SUBSTITUTION 

0INRLO3O1 


01 RB02 


1 .0 

IB 

01 RB01 


1.5 

IB 

01 RBOO 


3.0 

IB 

READ TABLE 

FOR SUBSTITUTION 

0PERTO5 


02 BTDD01 


1.0 


1.0 





TABLE A1.- CONTINUED 


02 BTD01 
1.0 

BBD 

1.0 

DBD 

02 BT01 
1.0 

BB 

1.0 

DB 

02 SINB01 
1.0 

BB 

1.0 

DB 

01 C0SB01 
1.0 

READ VARIABLES FOR COLLECTION OF COEFFICIENTS 

6MUCF12 
THO TH1 

TH2 BDDO BDD1 BDD2 BDO BD1 BD2 BO B1 B2 

READ TABLE 

FOR SUBSTITUTION 

@MULB04 
07 DBDD01 
1.0 

BDDO 

1.0 

BDD1 CSCY 

-2.0 

BD1 SNCY 

-1.0 

B1 CSCY 

1.0 

BDD2 SNCY 

2.0 

BD2 CSCY 

-1.0 

B2 SNCY 

05 DBD01 
1.0 

BDO 

1.0 

BD1 CSCY 

-1.0 

B1 SNCY 

1.0 

BD2 SNCY 

1.0 

B2 CSCY 

03 DB01 
1.0 

BO 

1.0 

B1 CSCY 

1.0 

B2 SNCY 

03 THTA01 
1.0 

THO 

1.0 

TH1 CSCY 

1.0 

TH2 SNCY 

READ TABLE 

FOR SUBSTITUTION 

@TRIG1201 
02 SNCY02 
0.5 
-.5 

C2CY 

01 SNCY01 CSCY01 

0.5 

S2CY 

02 CSCY02 

0.5 

0.5 

C2CY 

02 SNCY03 
0.75 

SNCY 

-.25 

S3CY 
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TABLE A 1 . - CONTINUED 


02 SNCY02 

CSCY01 

0.25 

CSCY 

-.25 

C3CY 

02 SNCY01 

CSCY02 

0.25 

SNCY 

0.25 

S3CY 

02 CSCY03 


0.75 

CSCY 

0.25 

C3CY 

03 SNCY04 


0.375 


-0.5 

C2CY 

0.125 

C4CY 

02 SNCY03 

CSDY01 

0.25 

S2CY 

-.125 

S4CY 

02 SNCY02 

CSCY02 

0.125 


-.125 

C4CY 

02 SNCY01 

CSCY03 

0.25 

S2CY 

0.125 

S4CY 

03 CSCY04 


0.375 


0.5 

C2CY 

0.125 

C4CY 


FORM MATRIX 
#TLAG#TFLP#TFLG 
FORM MATRIX 
#TFLG#POVE//PVEC 
MATRIX EXPRESSION 
03#PVEC 

% R 1 0 1 0 1 % R20201? R30301 

SUBSTITUTE TABLE INTO EXPRESSION 
@SUPR? R 1 % RX 

SUBSTITUTE TABLE INTO EXPRESSION 
@SUPR? R2? RY 

SUBSTITUTE TABLE INTO EXPRESSION 
@SUPR5 R3? RZ 
READ EXPRESSION 
?0MEZ01 
1.0 

DIFERENTIATE EXPRESSION 
TAU? RX? RXT 
DIFERENTIATE EXPRESSION 
TAU? RY? RYT 
DIFERENTIATE EXPRESSION 
TAU? RZ? RZT 
FORM EXPRESSION 
% RXD02 

1.0 % RXT 
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TABLE A1 . - CONTINUED 


-1.0 % RY^OMEZ 

FORM EXPRESSION 
% RYD02 

1.0 % RYT 

1.0 %omii RX 

FORM EXPRESSION 
% RZD01 

1.0 % RZT 

READ EXPRESSION 
/tSEBLOl 

0.5 KBT BT 

READ EXPRESSION 
JXVEL01 

1.0 MU CSCY 

READ EXPRESSION 
fYVELOl 

-1.0 MU SNCY 

READ EXPRESSION 
JZVEL01 

-1.0 LAMD 

FORM EXPRESSION 
HAVL102 

1.0 *XVEL 

-1.0 % RXD 

FORM EXPRESSION 
HAVL202 

1.0 %YVEL 

-1.0 t RYD 
FORM EXPRESSION 
*AVL302 

1.0 %ZVEL 

-1.0 % RZD 

READ MATRIX 
IVXYZ0301 
01 

1.0 *AVL1 

01 

1.0 %hVL2 

01 

1.0 *AVL3 

FORM MATRIX 
#TFLP#VXYZ#VX2Y C 
MATRIX EXPRESSION 
02IVX2Y 

% UT020 1 % U1P0301 

NEGATE EXPRESSION 
% UT 

NEGATE EXPRESSION 
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TABLE A1 CONTINUED 


% U1P 

SUBSTITUTE TABLE INTO EXPRESSION 
0DUMY* U1P? UP 
FORM EXPRESSION 
% MB02 

-0.5 % UT/& UT THTA RAC2 RB 

0.5 % UT/& UP RAC2 RB 

READ TABLE FOR SUBSTITUTION 
0DYNMO4O 1 
01 RAC201 RB01 
0.5 GAMI 

01 RAC201 RB02 
0.33333 GAMI 
01 RAC201 RB03 
0.25 GAMI 

01 RAC201 
1.0 GAMI 

DIFERENTIATE EXPRESSION 
BTD/J RXDySRXDQ 
DIFERENTIATE EXPRESSION 
BT* RXD/J RXQ 
DIFERENTIATE EXPRESSION 
BTD% RYD/tRYDQ 
DIFERENTIATE EXPRESSION 
BT? RYD£ RYQ 
DIFERENTIATE EXPRESSION 
BTDyt RZD/6RZDQ 
DIFERENTIATE EXPRESSION 
RZD/J RZQ 
READ EXPRESSION 
5tC0NS0 1 
1.0 

READ EXPRESSION 
£ CON 1 0 1 
- 1.0 


READ VARIABLES FOR COLLECTION OF COEFFICIENTS 
0PECFO6 

DBDD DBD DB BBDD BBD BB 
FORM LAGRANGIAN 
0803000100 
£SBEQ$PECF*E2D1 
% RXD? RXQ^CON 1 
% RYD* RYQ^CONI 
% RZD/6 RZQ/SCON 1 
% RXD/5RXDQ55CONS 
% RYD/tRYDQytCONS 
% RZD/tRZDQ/tCONS 
/tCONS/tSEBL/ICONS 
/JCONS/& MB^CONS 


@INRL§DUMY@PERT 
0INRL0DUMY0PERT 
0INRL0DUMY0PERT 
TAU0INRL0DUMY0PERT 
TAU@INRL§DUMY@PERT 
TAU@INRL@DUMY@PERT 
BT0DUM 1 @DUM 1 @PERT 


@DYNM§DUMY§PERT 
INITIALIZE MULTIBLADE 

03 MOPK CSCY SNCY C2CY S2CY C3CY S3CY C4CY S4CY C5CY S5CY C6CY S6CY C7CY S7CY 
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TABLE A1 . - CONCLUDED 


PERFORM MULTIBLADE TRANSFORMATION 
03 

@MULTB@TRIG@MUCF 
^SBEQJAEMI $AEM2/5AEM3 
END OF DATA 
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TABLE A2 . - SAMPLE OUTPUT 


READ VARIABLES FOR COLLECTION OF COEFFICIENTS 

RDCCVR DBDD DBD DB BBDD BBD BB 
FORM LAGRANGI AN 


a*##*#######*#**#*##*###***********##****# 


* 

* 

# 

* 


DETAILS OF THE EXPRESSION 
NUMBER OF TERMS 


» 

SBEQ * 

19 * 

* 


a##**####*###*#####*#####***********###### 


- GROUP OF TERMS WITH VARIABLE DBDD 


1 I.OOO* IB* 


- GROUP OF TERMS WITH VARIABLE DBD 


1 0.167*SNCY* MU*GAMI* 

2 0.125*GAMI* 


- GROUP OF TERMS WITH VARIABLE DB 


1 1.000* KBT* 

2 0.250*CSCY*SNCY* MU* MU*GAMI* 

3 0.167*CSCY* MU*GAMI* 

4 1.000* IB* 


- GROUP OF TERMS WITH VARIABLE BBDD 


1 1.000* IB* 


- GROUP OF TERMS WITH VARIABLE BBD 


1 0.167*SNCY* MU*GAMI* 

2 0 . 125*GAMI* 
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TABLE A2 . - CONTINUED 



- 

GROUP OF TERMS WITH VARIABLE 

BB 


1 

1.000* KBT* 



2 

0.167*CSCY* MU*GAMI* 



3 

1.000* IB* 



4 

0.250*CSCY*SNCY* MU* MU*GAMI* 


- 

REMAINING TERMS IN EQUATION 

- 

1 

-0.250*SNCY*SNCY*THTA* MU* 

MU*GAMI 

* 

2 

-0. 125*THTA*GAMI* 



3 

-0.333*SNCY*THTA* MU*GAMI* 



4 

0.250*SNCY* MU*LAMD*GAMI* 



5 

0. 167*LAMD*GAMI* 



PERFORM MULTIBLADE TRANSFORMATION 



MULTIBLADE EQUATION 1 



******************************************* 

* 



* 

* 

DETAILS OF THE EXPRESSION 

AEM1 

* 

* 

NUMBER OF TERMS 

17 

* 

* 



* 

******************************************* 


GROUP OF TERMS WITH VARIABLE 

THO 

- 

1 

-0.375* MU* MU*GAMI* 



2 

-0.375*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

TH1 

- 

1 

0.188*C3CY* MU* MU*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

TH2 

- 

1 

0.188*S3CY* MU* MU*GAMI* 



2 

-0.500* MU*GAMI* 
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TABLE A2 . - CONTINUED 


- GROUP OF TERMS WITH VARIABLE BDDO 


3.000* IB* 


- GROUP OF TERMS WITH VARIABLE BDO 


0.375*GAMI* 


- GROUP OF TERMS WITH VARIABLE BD2 


0.250* MU*GAMI* 


GROUP OF TERMS WITH VARIABLE BO 


1 3.000* IB* 

2 3.000* KBT* 



REMAINING TERMS IN EQUATION 


0.375* BBD*GAMI* 
3.000* BB* KBT* 
3.000* IB*BBDD* 
3.000* IB* BB* 
0 . 500*L AMD*GAMI * 





TABLE A2 . - CONTINUED 


MULTIBLADE EQUATION 2 



****************** ft*********************** 

* 



* 

* 

DETAILS OF THE EXPRESSION 

AEM2 

* 

* 

NUMBER OF TERMS 

19 

* 

* 



* 

****************************************** 

- 

GROUP OF TERMS WITH VARIABLE 


B 

1 0.188*C3CY* MU* MU*GAMI* 

- 

GROUP OF TERMS WITH VARIABLE 

TH1 

- 

1 

-0.094* MU* MU*GAMI* 



2 

-0. 188*GAMI* 



3 

-0.250*S3CY* MU*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

TH2 

- 

1 

0.250*C3CY* MU*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

BDD1 

- 

1 

1.500* IB* 



- 

GROUP OF TERMS WITH VARIABLE 

BD1 

- 

1 

0. 188*GAMI* 



2 

0.125*S3CY* MU*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

BD2 

- 

1 

3.000* IB* 



2 

-0.125*C3CY* MU*GAMI* 
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TABLE A2 . - CONTINUED 



- 

GROUP OF TERMS WITH VARIABLE 

BO 

- 

1 

0.250* MU*GAMI* 



2 

0.188*S3CY* MU* MU*GAMI* 



“ 

GROUP OF TERMS WITH VARIABLE 

B1 

- 

1 

1.500* KBT* 



2 

0.250*C3CY* MU*GAMI* 



“ 

GROUP OF TERMS WITH VARIABLE 

B2 

- 

1 

0.250*S3CY* MU*GAMI* 



2 

0. 188*GAMI* 



3 

0.094* MU* MU*GAMI* 



- 

REMAINING TERMS IN EQUATION 


- 

1 

0.188* BB*S3CY* MU* MU*GAMI* 


2 

0.250* BB* MU*GAMI* 



MULTIBLADE EQUATION 3 



****************************************** 

* 



K 

# 

DETAILS OF THE EXPRESSION 

AEM3 

# 

# 

NUMBER OF TERMS 

21 

* 

* 



* 

***** ************************************* 


GROUP OF TERMS WITH VARIABLE 

THO 

- 

1 

0.188*S3CY* MU* MU*GAMI* 



2 

-0.500* MU*GAMI* 



- 

GROUP OF TERMS WITH VARIABLE 

TH1 

- 

1 

0.250*C3CY* MU*GAMI* 
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TABLE A2 . - CONTINUED 


- GROUP OF TERMS WITH VARIABLE TH2 


1 -0.281* MU* MU*GAMI* 

2 -0. 188*GAMI* 

3 0.250*S3CY* MU*GAMI* 


- GROUP OF TERMS WITH VARIABLE BDD2 


1 

1 .500* 

IB* 




- 

GROUP OF 

TERMS 

WITH 

VARIABLE 


1 

0.250* 

MU*GAMI* 



- 

GROUP OF 

TERMS 

WITH 

VARIABLE 

BD1 

1 

-3.000* 

IB* 




2 

-0. 125*C3CY* 

MU*GAMI* 



GROUP OF 

TERMS 

WITH 

VARIABLE 

BD2 


1 0. 188*GAMI* 

2 -0.125*S3CY* MU*GAMI* 


- GROUP OF TERMS 

WITH 

VARIABLE 

BO - 

1 -0.188*C3CY* 

MU* 

MU*GAMI* 


- GROUP OF TERMS 

WITH 

VARIABLE 

B1 


1 0.250*S3CY* MU*GAMI* 

2 -0. 188*GAMI* 

3 0.094* MU* MU*GAMI* 
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TABLE A2 . - CONCLUDED 



- 

GROUP OF TERMS WITH VARIABLE B2 

1 

-0.250*C3CY* MU*GAMI* 

2 

1.500* KBT* 

- 

REMAINING TERMS IN EQUATION 

1 

0.250* BBD* MU*GAMI* 

2 

-0.188* BB*C3CY* MU* MU*GAMI* 

3 

0.375* MU*LAMD*GAMI* 

END OF 

DATA 
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TABLE A3.- SYMBOL DEFINITIONS 


Symbol 


FORTRAN Symbol 

Description 

A 


A 

B , C 

area of cross section, matrix name 
matrix names 

COS ip , cos 

2i> 

CSCY,C2CY 


cos 34>, cos 

4i|> 

C3CY,C4CY 


cos 6 , sin 

6 

COSB,SINB 


cos q, sin 

c 

COSZ , SINZ 


c d0 /a 


CDOA 

drag coefficient to lift curve slope 

E 


E 

Young's modulus 

*8 


IB 

flap moment of inertia 

K e 


KB 

flap spring stiffness 

r 


RB 

rotor radius (nondimensional) 

sin i|>, sin 

24* 

SNCY,S2CY 


sin 3f, sin 

4i|> 

S3CY,S3CY 


u T ,u p 


UT,UP 

tangential and perpendicular velocity components 
of flow over the airfoil 

u' 


US 

derivative of extensional displacement 

v' , w 1 


VS,WS 

bending slopes 

6 , 8,6 


BB,BBD,BBDD 

steady flap angles and derivatives 

68 , 66,66 


DB,DBD,DBDD 

perturbed flap angle and derivative 

8 o ,8 r B 2 


B0,B1 ,B2 

flap collective and cyclic angles 

B o ,B r B 2 


BD0,BD1 , BD2 

1st derivatives 

B o ,B r B 2 


BDD0.BDD1 , BDD2 

2nd derivatives 

u 


MU 

advance ratio 

0 


THTA 

pitch angle 

0^,0 ,0 
0 c' s 


TH0,TH1 ,TH2 
(THTO,THTC,THTS) 

collective and cyclic pitch angles 
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TABLE A3.- CONTINUED 


Symbol 

FORTRAN Symbol 

Description 

X 

LAMD 

uniform inflow 

q , n ,n 
x y z 

OMEX , OMEY , OMEZ 

expressions defining rotation vector 

*1 

GAMI 

Lock number 

pac/2 

RAC2 

air density 

* 


azimuth angle 

T 

TAU 

dimensionless time 

m b 

MB 

flap moment expression name 

c 

ZETA 

lag angle 

8,0,6 

BT , BTD , BTDD 

flap angle and time derivatives 

S ebl' S ebq 

SEBL ,SEBQ 

strain energy and single blade expression na 


PECF,MUCF 

group names for collection of coefficients in 
perturbed equations and multiblade equations 


MULB 

table containing multiblade relations 


TRIG 

table containing trigonometric relations 


PERT 

table containing perturbation relations 


DERV 

table containing differentiation rules 


INRL.DYNM 

tables containing integration rules 


SUPR 

table containing rules to suppress lag motion 


AEM1 ,AEM2,AEM3 

names of multiblade equations (expressions) 


TFLP,TLAG,POVE 

matrices giving transformation relations 

r ,r ,r 
x 7 y 7 z 

RX,RY,RZ 

expressions defining position vector 

r ,r ,r 
x 7 y 7 z 

RXT,RYT,RZT 

velocity expressions 


RXD , RYD , RZD 

a/3x {rxt,ryt,rzt} + Q x {rx,ry,rz} 


RXDQ,RYDQ,RZDQ 

3/3$ {RXD, RYD, RZD} 


RXQ,RYQ,RZQ 

3/36 {RXD, RYD, RZD} 
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TABLE A3.- CONCLUDED 


Symbol 

FORTRAN Symbol 

Description 


DUMY.DUM1 

tables containing dummy relations 


R1,LAEQ 

intermediate expression names 


E2D1 

name of the ordering scheme 


VX2Y 

intermediate matrix name 
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APPENDIX B 


QUASI -STATIC APPROXIMATION 


It is frequently possible to reduce the order of the system of equations 
describing the rotorcraft dynamics by considering a quasi-static approximation for 
certain degrees of freedom. Assume that the equations of motion have been reordered 
so that the quasi-static variables (Xq) appear last in the state vector. 


X 



( B 1 ) 


The quasi-static approximation consists of neglecting the acceleration and 
velocity terms of the quasi-static variables. Thus the equations of motion take the 
form 



(B2) 


The quasi-static variables are not described by differential equations but 
rather by linear algebraic equations. The solution for (Xq) then is simply 


= (a°°) 1 [-A® 1 X 1 - A^ 1 ^ - Aq 1 X 1 1 + B°V 


(B3) 


Substituting for Xq in the X^ equation of motion gives then the reduced 
order equations for the quasi-static approximation: 


[A* 1 - Aq°(A°°) _1 A° 1 ]X 1 + [a] 1 


c<c 


t«" - C<C>''O x i 


= [B 1 - Aq 0 ( Aq 0 ) _ ^B°]V (B4) 


In the present analysis, the quasi-static approximation is applied to the 
inflow dynamics and torsion degree of freedom of the blade. 
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TABLE 1.- COMMANDS AND BRIEF DESCRIPTION OF THEIR FUNCTION 


Command 

Function 


Input commands 

1. READ EXPRESSION 

Reads an expression name, number of terms, 
and term details. 

2. READ MATRIX 

Reads a matrix name, size, and term details 
for each element. 

3. READ TABLE FOR SUBSTITUTION 

Reads table name, number of relations, type 
of table, and relation details. 

4. READ DIFFERENTIATION TABLE 

Reads differentiation table name, number of 
relations, and relation details. 

5. READ GROUP AND ORDER OF 

Reads group of each variable and its order 

VARIABLES 

of magnitude. 

6. READ ORDERING SCHEME 

Reads order up to which variables of each 
group can be retained. 

7. READ VARIABLES FOR COLLECTION 

Reads a group name and variables for which 

OF COEFFICIENTS 

the coefficients have to be collected. 

8. READ TAPE 

Reads from a tape (disk), to separate 
steady, and perturbed terms. 

Algebraic manipulation commands 

9. FORM EXPRESSION 

Multiplies, adds, or subtracts two or more 
expressions. 

10. FORM MATRIX 

Multiplies two matrices including transpose 
multiplication. 

11. ADD MATRICES 

Adds or subtracts two compatible matrices. 

12. CONSTANT MATRIX 

Multiplies the matrix elements with a con- 
stant for expression having one term. 

13. MATRIX TO EXPRESSION 

Converts matrix elements to expression 
elements. 

14. SUBSTITUTE TABLE INTO EXPRESSION Substitutes a table of relations into an 

expression. 
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TABLE 1.- CONCLUDED 


Command 

Function 

Input commands - concluded 

15. DIFFERENTIATE EXPRESSION 

Differentiates an expression with respect to 
a variable. 

Application commands 

16. FORM LAGRANGIAN 

Uses Lagrangian formulation to generate 
governing equations. 

17. INITIALIZE MULTIBLADE 

Generates multiblade functions. 

18. PERFORM MULTIBLADE TRANSFORMATION Generates multiblade equations. 

19. WRITE TAPE 

Writes on tape (disk) compatible with sym- 
bolic program. 

20. NEGATE EXPRESSION 

Gives the negative of an expression. 

21. RESET COUNTER 

Erases unnecessary expressions. 

22. END OF DATA 

Stops execution. 

Program TRANS a 

Codes the equations into FORTRAN subrou- 
tines. 


a Executed separately at present. 








TABLE 2.- FORTRAN SYMBOLS 


Original 

variable 

FORTRAN 

symbol 

Original 

variable 

FORTRAN 

symbol 

R 

RAD 

v' 

VS 

u 

U 

v i 

LAMB 

u 

UD 

w 

W 

u T 

UT 

w' 

WS 

u p 

UP 

w 

WD 

Up 

UF 

X 

XCOR 

°F 

UFD 

0 

THTA 

V 

V 

4> 

PHI 

V 

VD 

G 

OMEG 
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TABLE 3.- TYPICAL INPUT TO HESL AND OUTPUT 
TO CALCULATE TANGENTIAL AND NORMAL 
VELOCITIES (U T AND U p ) 


READ MATRIX 

#LAFP0303 


03 



1.0 



-0.5 


VS VS 

-0.5 


WS WS 

01 



1.0 


VS 

01 



1 .0 


WS 

02 



-1.0 


VS^EONE 

-1 .0 


WS^ETWO 

03 



1.0 


%ETRE 

-1.0 


PHI THTA 

-1 .0 


VS WS THTA 

02 



1.0 


PHI 

1.0 


5&EFUR THTA 

02 



1.0 


VS/JETWO 

-1.0 


WS?EONE 

03 



-1.0 


2ETRE THTA 

-1.0 


PHI 

-1.0 


VS WS 

02 



-1.0 


PHI THTA 

1.0 


2EFUR 

READ MATRIX 

# VEL0301 


03 



-1.0 


UD 

1.0 


UFD 

1.0 


OMEG V 

02 



-1.0 


VD 

-1.0 


ztXUUF OMEG 

02 



-1.0 


OMEG RAD^LAMD 

-1.0 


WD 

FORM MATRIX 

#LAFP# 

VEL#AVEL*E2D1 

MATRIX 

EXPRESSION 

02#AVEL 


% UT0201* 

UP0301 

NEGATE 

EXPRESSION 

% UT 
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TABLE 3.- CONCLUDED 


NEGATE EXPRESSION 
% UP 


Note: #EONE, #ETWO, etc., are expressions 

read earlier in the program. 


Output of Uij. 


a#*##**#****#***##**#*##*#***#***#**#***##* 

* DETAILS OF THE EXPRESSION UT * 

* NUMBER OF TERMS 15 * 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


1 1.000* VS*0MEG* V* 

2 1.000* VD* 

3 1.000*XCOR*OMEG* 

4 1.000* U*OMEG* 

5 -1.000* UF*OMEG* 

6 -0.500* VS* VS* VD* 

7 -0.500* VS* VS*XC0R*0MEG* 

8 -0.500* PHI* PHI* VD* 

9 -0.500* PHI* PHI*XC0R*0MEG* 

10 -1.000* PHI*THTA* VD* 

11 -1.000* PHI*THTA*XCOR*OMEG* 

12 1.000* PHI*LAMB*OMEG* RAD* 

13 1.000* PHI* WD* 

14 1 .000*THTA*LAMB*0MEG* RAD* 

15 1.000*THTA* WD* 
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(a) Flow chart. 

Figure 1.- Aeroelastic analyses. 
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(b) Definition and storing of expressions and relations. 
Figure 1.- Concluded. 
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Figure 3.- Inflow dynamics. 
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HODGES AND DOWELL'S EQUATIONS, c = 5.0 

O COMPUTER GENERATED EQUATIONS, = 5.0 

□ HODGES AND DOWELL’S EQUATIONS, = 2.5 

O COMPUTER GENERATED EQUATIONS, = 2.5 



Figure 4.- Comparison of lead-lag damping from two sets of equations for varying 
pitch angle 0 for soft in-plane blade for two torsional frequencies, io v = 0.7, 
o, w = 1.15, R = 1.0, a = 0.1. 
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Figure 6.- Lead-lag mode damping versus pitch angle 0, from flap-lag analysis for a 
rigid blade model and an elastic blade model for a soft in-plane blade, w v = 0.7, 
u u = 1.15, R = 1.0, a = 0.05, C d = 0.01. 



Figure 7.- Lead-lag mode damping versus pitch angle 9, from flap-lag analysis for 
rigid blade model and an elastic blade model for a stiff in-plane blade, 


0 ). 


= 1.5, <» w = 1.15, R = 1.0, o = 0.05, C ri = 0.01 



-30 r 



Figure 8.- Lead-lag regressing mode damping with dynamic inflow versus pitch angle 
0, from rigid blade model and elastic blade model for a soft in-plane blade, 
u> v = 0.7, u> w = 1.15, R = 1.0, a = 0.05. 
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Figure 9.- The effect of torsion frequency and dynamic inflow on lead-lag regressing 
mode damping versus lead-lag frequency, 9 = 0.3 rad, R = 0.0, tu w = 1.15, a = 0.1. 
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Figure 10.- The effect of torsion frequency and dynamic inflow on lead-lag regress- 
ing mode damping versus lead-lag frequency, 0 = 0.3 rad, R = 1.0, u) w = 1.15, 
o=0.1. 


87 




Figure 11.- The effect of elastic coupling and dynamic inflow on lead-lag regressing 
mode damping versus lead-lag frequency, 9 = 0.3 rad, = 5.0, io w = 1.15, 

= 0 . 1 . 


WITH DYNAMIC INFLOW 

WITHOUT DYNAMIC INFLOW 



0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0 
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Figure 12.- The effect of torsion frequency and dynamic inflow on lead-lag regress 
ing mode damping versus elastic coupling, 8 = 0.3 rad, u > v = 1.4, u> w = 1.15, 

= 0 . 1 . 


WITH DYNAMIC INFLOW 

WITHOUT DYNAMIC INFLOW 



Figure 13.- The effect 
mode damping versus 


of torsion frequency and dynamic inflow on flap regressing 
elastic coupling, 0 =0.3 rad, a> v = 1.4, u> w = 1.15, a = 0.1. 
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Figure 14.- Dynamic inflow modeling effects on lead-lag regressing mode damping 
versus elastic coupling. e = 0.3 rad, = 5.0, u> v = 1.4, w w = 1.15, o = 0.1. 
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WITH DYNAMIC INFLOW 


WITHOUT DYNAMIC INFLOW 



Figure 15.- Lead-lag regressing mode damping from flap-lag-torsion and quasi-static 
torsion model with and without dynamic inflow versus elastic coupling parameter, 
<*>,), = 5.0, w w = 1.15, u> v = 1.4, o = 0.1, 0 = 0.3 rad . 
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Figure 16.- Lead-lag regressing mode damping from a quasi-static torsion model and a 
rigid blade flap-lag model with kinematic pitch-flap and pitch-lag couplings, 

0 = 0.3 rad , w v = 1.4, oi w = 1.15, = 5.0. 




C - 
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AJ2R = Vsin a + Vj X 

SIDE VIEW TOP VIEW 


Figure 17.- Aerodynamic velocity components in forward flight. 
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Figure 18.- The effect of number of degrees of freedom used in trim analysis on 
lead-lag damping versus advance ratio, soft in-plane, w v = 0.7, u w = 1.15, 

= 3.0, a = 0.1, R = 1.0, propulsive trim. 
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ADVANCE RATIO, fi 


Figure 19.- The effect of number of degrees of freedom used in trim analysis on 
lead-lag damping versus advance ratio, stiff in-plane, u = 1.4, u = 1.15, 
m. = 3.0, a = 0.1, R = 1.0-, propulsive trim. 




Figure 20.- Lead-lag mode damping versus advance ratio for a flap-lag-torsion model 

and a flap-lag model, soft in-plane, u> = 0.7, = 1.15, u>. = 3.0, a = 0.1, 

R = 1.0. v w <j, 
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NUMBER OF MODES 



Figure 21.- Lead-lag mode damping versus advance ratio for a flap-lag-torsion model 
and a flap-lag model, soft in-plane, « = 0.7, ui u = 1.15, a>. = 3.0, a = 0.1, 

R = 0.0. v w $ 
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Figure 22.- Lead-lag mode damping versus advance ratio for a flap-lag-torsion model 
and a slap-lag model, stiff in-plane, u> v = 1.4, <u w = 1.15, = 3.0, a = 0.1, 

R = 1 .0. 


o 


NUMBER OF MODES 
2 FLAP, 2 LAG, 2 TORSION 


a 2 FLAP, 2 LAG, 1 TORSION 



Figure 23.- Lead-lag mode damping versus advance ratio for a flap-lag-torsion model 

and a flap-lag model, stiff in-plane, w = 1.H, u = 1.15, w. = 3.0, a = 0.1, 
p = n n v w <p 




Figure 24.- Lead-lag mode damping versus elastic coupling for a flap-lag-torsion 
model and flap-lag model, stiff in-plane, u = 1.4, w = 1.15, w. = 3.0, a = 0.1 
y = 0.25. 


1 01 


} 



Figure 25.- The effect of torsion and dynamic inflow on lead-lag regressing mode 
damping versus advance ratio, soft in-plane, co v = 0.7, u w = 1.15, = 3.0, 

o = 0.1, R = 1.0. 
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Figure 26.- The effect of torsion and dynamic inflow on lead-lag regressing mode 
damping versus advance ratio, soft in-plane, u = 0.7, o> = 1.15, w, = 3.0, 
o = 0.1, R = 0.0. 9 



O FLAP-LAG-TORSION 



Figure 27.- The effect of torsion and dynamic inflow on lead-lag regressing mode 
damping versus advance ratio, stiff in-plane, w v = 1.4, u> w = 1.15, = 3.0, 

a = 0,1, R = 1,0. 




Figure 28.- The effect of torsion and dynamic inflow on lead-lag regressing mode 


damping versus advance ratio, stiff in-plane, 
a = 0.1, R = 0.0. 


1.4, u M = 1.15, 


3.0, 




0 .2 .4 .6 .8 1.0 
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Figure 29.- The effect of torsion and dynamic inflow on lead-lag regressing mode 
damping versus elastic coupling, stiff in-plane, w v = 1.4, u» w = 1.15, = 3.0, 

o=0.1, y=0. 25. 
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WITHOUT DYNAMIC INFLOW 



ADVANCE RATIO, n 


Figure 30.- Comparison of constant coefficient approximation and Floquet analysis, 
stiff in-plane, a> v = 1.4, u> w s 1.15, u = 3.0, o = 0.1, R = 1.0. 
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